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We present a numerical code for calculating the local gravitational self-force acting on a pointlike 
particle in a generic (bound) geodesic orbit around a Schwarzschild black hole. The calculation is 
carried out in the Lorenz gauge: For a given geodesic orbit, we decompose the Lorenz-gauge metric 
perturbation equations (sourced by the delta-function particle) into tensorial harmonics, and solve 
for each harmonic using numerical evolution in the time domain (in l-fl dimensions). The physical 
self-force along the orbit is then obtained via mode-sum regularization. The total self-force contains 
a dissipative piece as well as a conservative piece, and we describe a simple method for disentangling 
these two pieces in a time-domain framework. The dissipative component is responsible for the loss 
^ ' of orbital energy and angular momentum through gravitational radiation; as a test of our code we 

\ demonstrate that the work done by the dissipative component of the computed force is precisely 

balanced by the asymptotic fluxes of energy and angular momentum, which we extract independently 
from the wave-zone numerical solutions. The conservative piece of the self-force does not affect the 
time-averaged rate of energy and angular-momentum loss, but it influences the evolution of the 
orbital phases; this piece is calculated here for the first time in eccentric strong-field orbits. As a 
\^ • first concrete application of our code we recently reported the value of the shift in the location and 

frequency of the innermost stable circular orbit due to the conservative self-force [Phys. Rev. Lett. 
102, 191101 (2009)]. Here we provide full details of this analysis, and discuss future applications. 
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The prospects for detecting gravitational waves from the inspiral of compact objects into massive black holes have 
motivated, over the past decade, research in effort to understand the general-relativistic orbital evolution in such 
systems. The underlying elementary theoretical problem is that of a pointlike mass particle in a strong-field orbit 
around a Kerr black hole of a much larger mass. The dynamics of such systems can be described in a perturbative 
^■f^ . fashion in terms of an effective gravitational self- force (SF) I'-'s'l; knowledge of this force is a prerequisite for describing 
^sj ' the precise evolution of the orbit and the phasing of the emitted gravitational waves. There is an active research 
\ program focused on the development of computational methods and actual working codes for the SF in Kerr spacetime 
■ [6|. This research agenda is being pursued in incremental steps, through exploration of a set of simplified model 
' problems with increasing complexity and physical relevance. Much of the initial work has concentrated on a scalar- 
^ \ field toy model 0-0, but more recently workers have begun to tackle the gravitational case The state 

't^ ' of the art is represented by three independent calculations of the gravitational SF for circular geodesic orbits in 
• ^ Schwarzschild geometry [3- These calculations use different analytic and numerical methods (and they even 

invoke different physical interpretations of the SF) , but they were shown to be fully consistent with each other [l8, To] . 
J-j These calculations were also shown to be consistent with results from post-Newtonian theory in the weak-field limit 



In the current work we extend the analysis of Ref. flE\ (hereafter "Paper I") from the special class of circular 
geodesies to generic (bound) geodesies of the Schwarzschild geometry. This generalization is astrophysically relevant 
because real inspirals often remain quite eccentric up until the eventual plunge into the massive hole |22| . At a 
more fundamental level, the generalization to eccentric orbits is interesting because it allows us to start exploring 
in earnest the conservative effects of the SF — for instance, how it influences the orbital precession. Eccentric orbits 
have already been considered in calculations of the scalar [l^] and electromagnetic (EM) [23| SFs by Haas. While 
these calculations are of a less direct astrophysical relevance, they offer an important test-ground for computational 
techniques potentially applicable in the gravitational problem too. Indeed, many elements of our numerical method 
take their inspiration from Haas' work. 

The numerical code we present here takes as input the two orbital parameters of an eccentric Schwarzschild geodesic 
(the semi-latus rectum and eccentricity, to be defined below) , and returns the value of the Lorenz-gauge gravitational 
SF along this geodesic. The dissipative and conservative pieces of the SF are returned separately. Here we do not 
consider the evolution of the orbit under the effect of the SF, but leave this important next step for future work. We 
envisage using, to this end, a version of the "osculating geodesies" method [l^], which takes as input the value of the 
SF along geodesies tangent to the actual inspiral orbit. A systematic framework for analyzing the long-term evolution 
of the inspiral orbits, using multiple-scale perturbation methods, was recently developed by Hinderer and Flanagan 
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[H (cf. Sec. VII of Gralla and Wald H). 

Our strategy is sirailar to that of Paper I. Its basic elements are (i) the Lorenz-gauge perturbation formahsm of 
Barack and Lousto [2^ , (ii) a finite-difference algorithm for numerical integration of the Lorenz-gauge perturbation 
equations in the time domain, and (iii) mode-sum regularization [27l - [30| . The perturbation formalism is based on a 
tensor-harmonic decomposition of the perturbed Einstein equations in the Lorenz gauge. The equations are augmented 
with "gauge damping" terms designed to suppress gauge violations [1^, and are written as a set of 10 hyperbolic 
equations (for certain linear combinations of metric components) which do not couple at their principal parts. These 
equations are sourced by the (tensor-harmonic modes of the) particle's energy-momentum, modeled with a delta- 
function distribution along the specified eccentric geodesic. The equations are solved numerically mode by mode in 
the time domain using characteristic coordinates on a uniform l-|-l-dimensions mesh. The non-radiative monopole 
and dipole modes cannot be evolved stably in this manner; instead, we solve for these two modes separately in 
the frequency domain, using the recently introduced "extended homogeneous solutions" technique [3l| to cure the 
irregularity of the Fourier sum near the particle. The code records the value of the perturbation modes and their 
derivatives along the orbit (each mode has a C° behavior at the particle and hence a well-defined value there, as well 
as a well-defined "one-sided" derivatives) . These values are then fed into the "mode-sum formula" [295 : which returns 
the physical SF through mode-by-mode regularization. 

One of the primary advantages of the time-domain approach is that eccentric orbits — even ones with large 
eccentricity — are essentially "as easy" to deal with as circular orbits, with computational cost being only a weak 
function of the eccentricity [32| . Also, a time-domain code for circular orbits can be upgraded with relative ease to 
accommodate eccentric orbits (such a generalization is radically less straightforward in the frequency domain) . Still, 
there are several important technical issues which arise in the time-domain upgrade from circular to eccentric orbits, 
and need to be addressed. We list some of these issues below. 

• Most obvious, the computational burden increases significantly because the parameter space for geodesies turns 
from ID (circular) to 2D (eccentric). Moreover, for each given geodesic parameters the SF becomes a function 
along the orbit (it has a constant value along a circular geodesic), and one is required to obtain this function 
over an entire radial period. The latter becomes a technical hurdle in situations where the radial period is very 
large — e.g., close to the last stable orbit, or for orbits with very large radii. 

• In paper I we were able to improve the convergence rate of our finite-difference algorithm using a Richardson-type 
extrapolation to the limit of a vanishing numerical grid-cell size. That was possible because in the circular-orbit 
case the numerical mesh could be easily arranged such that the local discretization error varied smoothly along 
the orbit. This cannot be achieved in any simple way when the orbit is eccentric, and as a result one cannot 
implement a similar Richardson extrapolation. The practical upshot is that one is forced to implement a higher- 
order finite-difference scheme: a 2nd-order-convergent algorithm (as in Paper I) proves insufficient in practice. 
For this work we developed an algorithm with a 4th-order global convergence. The algorithm takes a rather 
complicated form near the particle's trajectory, where the field (the Lorenz-gauge metric perturbation) has 
discontinuous derivatives. To somewhat lessen this complexity (and reduce the number of grid points needed as 
input for the finite-difference formula) the algorithm makes use of suitable junction conditions across the orbit. 
The eventual numerical scheme is considerably more sophisticated — and involved — compared to that of Paper 
I. 

• In the mode-sum scheme one first calculates the contribution to the "full" (pre-regularization) force from each 
tensorial-harmonic mode of the perturbation, and then decomposes this into spherical harmonics. The necessary 
input data for the mode-sum formula are the individual spherical-harmonic contributions. This procedure 
involves the implementation of a tensor-scalar coupling formula, whose details depend on the orbit in question. 
The coupling formula simplifies considerably in the circular-orbit case; it reverts to its full complicated form 
[Eq. (|2.25p with Appendix [C] when eccentric orbits are considered. 

• The computation of the monopole and dipole contributions to the SF (which we perform in the frequency 
domain, as mentioned above) becomes much more involved in the eccentric-orbit case. First, the spectrum 
of the orbital motion now includes all harmonics of the radial frequency, and one has to calculate and add 
up sufficiently many of these harmonics. A second, more technically challenging complication arises from the 
fact that the perturbation becomes a non-smooth function of time across the orbit (at a given radius), which 
disrupts the high-frequency convergence of the Fourier sum at the particle (a behavior reminiscent of the Gibbs 
phenomenon). A general method for circumventing this problem in frequency-domain calculations was devised 
recently in Ref. [3l|, and we implement it here for the first time. 

• In exploring the physical consequences of the SF it is useful to split the SF into its dissipative and conservative 
pieces, and discuss their corresponding effects in separate. This splitting is straightforward in the circular-orbit 
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case: The conservative piece is precisely the (Schwarzschild) r component of the SF, while the (Schwarzschild) 
t, (p components exactly account for the entire dissipative effect. This is no longer true for eccentric orbits, 
where each of the Schwarzschild components mixes up both dissipative and conservative pieces, and it is not 
immediately obvious how to extract these pieces individually. Here we suggest and implement a simple new 
method for constructing the dissipative and conservative pieces out of the computed Schwarzschild components 
of the SF (without resorting to a calculation of the advanced perturbation). The method takes advantage of 
the general symmetries of Schwarzschild geodesies. 

With the computational framework in place, we can start to explore the physical effects of the gravitational SF. In 
this article we concentrate on two such effects. First, we calculate the loss of orbital energy and angular momentum, 
over one radial period, due to the dissipative piece of the SF. We extract these quantities directly from the computed 
SF along the geodesic orbit (for a sample of orbital parameters) . These "lost" energy and angular momentum must be 
balanced by the total amount of energy and angular momentum in the gravitational waves radiated to spatial infinity 
and into the black hole over a radial period. We derive formulas for extracting these quantities from the far-zone 
and near-horizon numerical Lorenz-gauge solutions, and demonstrate numerically that they agree well with the values 
computed from the local SF. Our values for the energy and angular momentum losses also agree with those previously 
obtained by others using other methods. 

The second effect we consider is conservative, and cannot be inferred indirectly from the asymptotic gravitational 
waves: It is the conservative shift in the location and frequency of the Innermost Stable Circular Orbit (ISCO). The 
analysis of the ISCO shift requires knowledge of the SF along slightly eccentric geodesies near the last stable orbit, and 
our code provides the necessary SF data for the first time. We reported the results in a recent Letter [s^, and here we 
describe our analysis in full detail. The quantitative determination of the ISCO shift is important in that it provides 
a strong-field benchmark for calibration of approximate (e.g., post- Newtonian) descriptions of binary inspirals. Our 
result for the ISCO frequency shift has already been incorporated by Lousto et al. in their "empirical" fitting formula 
for predicting the remnant mass and spin parameters in binary mergers (s^ . [35| ; and by Damour p]| for breaking the 
degeneracy between certain unknown parameters of the Effective One Body (FOB) formalism. 

Perhaps of a more direct relevance to the problem of the phase evolution in binaries with extreme mass-ratio is the 
effect of the SF on the periapsis precession of the eccentric orbit — also a conservative effect. SF corrections to the 
precession rate have been analyzed for weak- field orbits and within the toy model of the EM SF [s^ [s^l , but never 
before for the gravitational problem in strong field. Our code generates the SF data necessary to tackle this problem 
for the first time. We leave the detailed analysis of SF precession effects to a forthcoming paper. 

The paper is organized as follows. In Sec. II we review the relevant theoretical background: bound geodesies in 
Schwarzschild geometry, the Lorenz-gauge metric perturbation formulation, and the construction of the SF via the 
mode-sum formula. Section III describes our numerical method in detail, and in Sec. IV we present numerical results 
for a few sample eccentric orbits, including a "zoom-whirl" orbit. We explain how the dissipative and conservative 
pieces of the computed SF can be extracted from the numerical data, and present these two pieces separately in a few 
sample cases. We also analyze the dissipative effect of the SF and demonstrate the consistency between the dissipated 
energy and angular momentum inferred from the local SF, and that extracted from the asymptotic gravitational 
waves. Section V covers the ISCO-shift analysis, and in Sec. VI we summarize and discuss future applications of our 
code. 

Throughout this work we use standard geometrized units (with c = G = 1), metric signature — I — I — h, and (unless 
indicated otherwise) Schwarzschild coordinates = {t,r,0,ip). 



II. THEORETICAL BACKGROUND 
A. Eccentric geodesies in Schwarzschild geometry 



In this work we consider a pointlike particle with mass /i in a bound orbit around a Schwarzschild black hole of 
mass M ^ /i. In the limit /i — > the trajectory of the particle is a timelike geodesic of the background Schwarzschild 
spacetime. We parameterize this geodesic by proper time r, in the form x!^{t) — (tp(T), rp(T), 6'p(t), (/5p(r)), with 
corresponding four-velocity u'' = dz'^/dr. Without loss of generality we take 0p(r) = 7r/2. The geodesic equations of 
the particle are given in integrated form as 



(2.1) 



dt-p £ dipp C 

dr f(rp)' dr r^' 

Rirp,C% R{r,C') ^ f{r) (l + ^ ) , (2.2) 
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where / = 1 — 2M/r, and £ = —Ut and C = are the integrals of motion corresponding to the particle's specific 
energy and angular momentum. 

When > 12Af^ the effective potential of the radial motion, R{r,C^), has a maximum and a minimum and 
hence eccentric (bound) orbits exist. These orbits can be parametrized by the two values of Vp at the turning points, 
Tmin and Tinax ( "pcriastron" and "apastron", respectively). We may alternatively parameterize the orbits by the 
(dimensionless) semi-latus rectum p and eccentricity e, defined through 



2r, 



P 



mill' max 



M(r„ 



(2.3) 



From the two conditions R{r„ 



Ri^max) ~ one readily obtains 



(p-2-2e)(p-2 + 2e) 
p{p — 3 — e^) 



(2.4) 



Bound geodesies have < e < f and p > 6 + 2e 38]. Points along the separatrix p = 6 + 2e (where £^ equals the 
maximum of the effective potential) represent marginally unstable orbits. Stable circular orbits are those with e = 
and J3 > 6, for which £^ equals the minimum of the effective potential. The point (p, e) = (6, 0) in the e-p plane, 
where the separatrix intersects the e = axis, is referred to as the innermost stable circular orbit (ISCO) — see Fig. 

m 
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FIG. 1: Parameter space for bound geodesies in Schwarzschild spacetime. The (dimensionless) "semi-latus rectum" p and 
"eccentricity" e are defined in Eq. (|2.3|l . Bound geodesies have e > and p > 6 + 2e. Points along the separatrix p = 6 + 2e 
represent marginally unstable orbits. Stable cireular orbits lie along the axis e = for p > 6. The point (p, e) = (6, 0) is the 
ISCO. 



Following Ref. [3g, we introduce the monotonically increasing "radial phase" parameter x, defined so that the 
radial motion obeys 



^p(x) 



pM 



f + e cos X 



(2.5) 



Note X = 27rn (n integer) correspond to periastron passages. In terms of x^ the t and Lp components of the geodesic 
equations (|2.ip are re-expressed as 



dip _ Mp^ 

dx (p — 2 — 2ecosx)(l + ecosx)^ 



dx 



P 

p — 6 — 2ecosx' 



' (p-2^2e)(p-2-f2e) 
p — 6 — 2ecosx ' 



(2.6) 
(2.7) 



and the radial velocity reads 



— £esmx\ 



' p — 6 — 2e cos X 
{p-2-2e){p-2 + 2e)' 



(2.8) 
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The period of the radial motion can be derived by integrating Eq. (I2.6P with respect to x- 

T^^ ^dx. (2.9) 
Jo ax 

With the initial conditions — ip^ — a.t x — the particle's geodesic trajectory is fully specified by Eqs. (|2.5p . 
p.6p and (|2.7p . The functions tp(x) and </'p(x) cannot be written explicitly in analytic form, but it is easy to obtain 
them numerically, for given p and e, at any desired accuracy. 



B. Gravitational self- force via mode-sum regularization 

When ji is finite (yet still much smaller than M), the particle experiences a gravitational SF, F°'[^ 0{^^)], and the 
equation of motion is formally given by 

Here we use ip(f) to denote the (non-geodesic) trajectory under the effect of the SF, with f representing proper time 
along this trajectory and vf^ = dx^/df. The covariant derivatives D/Df are taken with respect to the background 
geometry. From symmetry we have = 0. Furthermore, assuming the normalization Uau" = —1, we have the 
orthogonality condition UaF°' = 0, which interrelates the remaining 3 components of the SF. All in all, then, there 
are two non-trivial independent components of the SF to be determined. 

In this work we do not consider the evolution of the orbit under the effect of the SF, i.e., we do not seek to obtain 
consistent solutions of Eq. (|2.10p . Rather, we are interested in calculating the value of the SF F°'{t) along a fixed, 
geodesic orbit x'^{t), with given values of p, e. We envisage that the SF information F"(t;p, e) (calculated over the 
space of p, e) could be used, in a second step, to calculate the orbital evolution in situations where at any given time 
the orbit deviates only very slightly from a geodesic of the background, and the evolution takes place over a timescale 
much longer than the radial period ("adiabatic approximation"). Here, however, we concentrate on the first step, 
leaving the investigation of the orbital evolution for future work. 

The gravitational SF acting on the particle at any a point along the geodesic x'^{T]p,e) is calculated using the 
mode-sum formula l27H29ll 



= E [^fun± -AIL-B'^]^Y. (2-11) 



1=0 1=0 

Here L = I + 1/2, and Ff^'j are the multipole modes of the "full" force field constructed from the (physical, retarded) 
Lorenz-gauge metric perturbation as prescribed in Sec. IHD] below. The subscript ± refers to the two possible values 
of J^fu'ii at 2;p, resulting from taking one-sided radial derivatives of the metric perturbation from either r — > r+ or 
r ^> Tp . and B" are the "regularization parameters" , given by [1^ [s^ 

^± = T^, ^± = T^, Al^O, (2.12) 



= 

B'' = 

Bf = 



7rr2[/3/2 



-K{w) +2{1 - U)E{w) 



(f 2 + /pC/) k{w) - [2£2(i -u)- /pC/(l - 2U)\ E{w) 



K{w)- (1 + 2^ ]E{w) 



(2.13) 



where K{w) = /J^^^(l — w sin'^ x) ^^"^dx and £^(1;;) = /J^^^(l — w sin^ x^l'^dx are complete elliptic integrals of the first 
and second kind, respectively, fp = f{rp) = 1 — 2M/rp, and 



£2 



w 



U = 1 



£2 



(2.14) 
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It is important to remember that the gravitational SF [as also the trajectory Xp{T) itself] is a gauge-dependent entity 
[3^ . The mode-sum formula (|2.1ip is formulated in the Lorenz-gauge, and requires as input the modes F^^^ derived 
from the Lorenz-gauge metric perturbation. In our approach this perturbation is obtained by tackling the linearized 
Lorenz-gauge Einstein equations directly, in the time domain. We proceed by reviewing the relevant Lorenz-gauge 
perturbation formalism. 



C. Metric perturbation in Lorenz gauge 

Ref. (26| presented a formulation of the Lorenz-gauge perturbation equations in Schwarzschild spacetime, amenable 
to numerical treatment in the time domain. In Paper I we applied this formulation (with some minor modifications) in 
our study of circular orbits. Here we shall use the same Lorenz-gauge formulation to obtain our metric perturbation, 
and we describe it here as applied to generic eccentric orbits. 

Let Qfj^i, be the metric of the background Schwarzschild geometry, and /i^^ be the physical (retarded) metric per- 
turbation due to the particle moving on the geodesic Xp (r;p, e). We assume hfj,i, is given in the Lorenz gauge, i.e., it 
satisfies 

V'"^ = - 0, (2.15) 

where hf^^, = hf^^, — (l/2)g^^h is the trace-reversed perturbation. The corresponding Einstein equations, linearized in 
h^i, over g^^, take the compact form 

+ 2R'^/,ha,0 = -167rr^., (2.16) 

where a semicolon denotes covariant differentiation with respect to 5^11/, and h = g^'^h^^, is the trace of /i^^. T^j^ on 
the right-hand side is the energy-momentum tensor associated with the particle, given by the distribution 

T^^ = /i / = i— dr, (2.17) 



where g is the determinant of (7^1^. 

It is well known that the hyperbolic set (|2.16p admits a well-posed initial-value formulation, and that the gauge 
conditions ()2.15|) are satisfied automatically if only they are satisfied on the initial (Cauchy or characteristic) surface. 
However, in a time-domain numerical implementation of Eqs. (j2.16l) it is usually impossible to satisfy the gauge 
conditions (j2.15p precisely on the initial surface, and, even if one succeeded to do so, finite differencing errors would 
usually lead to an uncontrollable violation of the gauge conditions (unless one somehow actively imposes the gauge 
conditions during the evolution). This problem, and its resolution, were discussed in Ref. 26], and we follow here the 
same method. To the original field equations (|2.16p we add "divergence dissipation" terms, in the form 

+ 2i?"^^^/i„^ + f{t,X + t^Z^) = -WnT^,, (2.18) 

where = (1, /^^, 0, 0) and ~ {fZr, Zr, Zg, Z^). While the inclusion of these extra terms does not (in principle) 
affect the solutions of the equations, it guarantees that violations of the gauge condition arc efficiently damped during 
the numerical time evolution. [This can be shown by considerin g th e divergence of Eq. (|2.18p . noticing that this yields 
a hyperbolic equation for Z^, with a manifest dissipation term |26|.] 

Owing to the spherical symmetry of the background geometry, the field equations (j2.18p are separable into tensorial 
spherical-harmonics using the substitution 

00 I 10 

V = E E«F'^^^^""(^i)^i^^''"(^,^;0 (2.19) 

1=0 m=-l 1=1 

(and similarly for the source T^^). The tensorial-harmonic basis vlf'J''™' and normalization factors (i = 1, . . . , 10) 

are the ones defined in Ref. [1^, except that here (as also in Paper I) we take F^^?'™ — )■ /(r) x V^^?'™ (this modification 
is needed to achieve h^'^^'-^ — ^const as r — 2A/, in line with the behavior of all other functions /i^*)'"*). In Appendix 
|A]we give explicit formulae for reconstruction of the various Schwarzschild components h^i, out of the 10 scalar-like 
functions W^'^"'{r,t). 

The above substitution reduces the field equations (|2.18p to the coupled set of 2-dimensional hyperbolic equations 

□;iW'™+7\/(«J/iO)'™ ^5'W'™^(r_rp) (i = 1,...,10). (2.20) 
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Here a box represents the 2-diniensional scalar-field wave operator 

/ 



U = duv + V{r), V{r) 



4^2 



2M 



(2.21) 



where v and u are the standard Eddington-Finkelstein null coordinates, defined through v = t + r^, and u = t — r^,, 
with = r + 2Mln[r/(2M) - 1]. The terms A^^(']j/i(-')'" (summation over j implied) involve first derivatives of 

the /i(^''™'s at most — hence the principal part of the set (I2.20p is entirely contained in the term D/it')'™. are 
the source terms for the point particle, constructed from the tensor- harmonic coefficients of T^^. In Appendix IB] we 
give explicit expressions for both A^*-^*^.'/i(-')'™ and S**-*^'™. The time-radial functions /i^*)'™ also satisfy four elliptic 
equations, which arise from the Lorenz gauge conditions ()2.15p . These relations, too, are given in Appendix iBl 



D. Construction of the full force mode 



Given the Lorenz-gauge metric perturbation hap, the full force modes ffyn^ appearing in the mode-sum formula 
(|2.1ip are formally constructed as we now prescribe. 

First, following we define the "full- force field" as a tensor field at arbitrary spacetime points a;, for a given 
(fixed) worldline point Xp (where the SF is to be calculated): 

F,Z,,{x; Xp) = fik^P^'ix; Xp)hp,,s- (2.22) 

Here the trace-reversed metric perturbation, /iq/j, is evaluated at x, and 

k'^'^''\x; Xp) = g'^^u^u''/2 - g'^^^u^'u^ - u"u''uV/2 + u"/''uV4 + 5" V/4, (2.23) 

where g°'^ is the background metric at x, and are the values of the contravariant Schwarzschild components of the 
four- velocity at Xp (treated as fixed coefhcients) . In principle, one can choose to extend the quantity fc^'^T* off the 
worldline in any one of many different ways (a few natural choices are discussed in Ref. [2§|). The specific choice made 
here is advantageous in that it guarantees a finite mode-coupling in Eq. (|2.25p below. Our choice of extension does not 
correspond to any of the choices made in ^29,], but it can be shown (using the methods of 29]) that the rcgularization 
parameters associated with our extension are the same as those of the "fixed contravariant components" extension 
defined in [2^ — these are the parameters whose values we state above in Eqs. (|2.12p and (|2.13p . 

In the next step we expand hap in tensor harmonics as in Eq. (|2.19p and substitute in Eq. (|2.22p . Taking the limits 
r — Tp and t — ^ tp (but maintaining the 0, (f dependence), the full force takes the form 

2 oo I 

[FU{6,^-rp,tp)]^ = ^ {/ot">^'™ + /fi"sin2(?r'™ + /2"i™cosf?sin0y^'" 

P 1=0 m=-l 

+ f^i^ sin^ e y'™ + /4 (cos eiy'" - sin ev^^) 

+ sin e Y^^ + ff^ sin^ 9 Y^^ + f^^ cos 9 sin^ 9 F'J^ } , (2.24) 

where Y^"^{9,tf) are the spherical harmonics, and the coefficients /^j.™ ^'"^ constructed from the perturbation fields 
J^(i)im their first r and t derivatives, all evaluated at Xp. The labels +/— correspond to taking one-sided derivatives 
from r — >■ or Tp , respectively. The explicit expressions for the 's are shown in Appendix [C] 

Since the mode-sum formula (12. lip requires as input the spherical harmonic modes of the full force, we must now 
re-expand Eq. (|2.24p in terms of spherical harmonics. With the help of the identities given in Appendix iDl we obtain 

2 ' 

jpal _ M \ ^ vlm/n ,„ ^ „ / Tr"'^3'™ i -r-al-2,m . -j^ai-l,m , -r-ahn , -r-al + ljn . -j^ai+2,m . -^al+S.mX ro OK\ 
-^fulli - ^ 2^ ^ It'pi'PpJX |-^(-3) +-^(-2) +-^(-1) +-^(0) +-^( + 1) +-^(+2) +-^(+3) /_^' ('^■^^) 

P m= — l 

where each of the functions is a certain linear combination of the /"±"'s (with the same l,m) — the explicit 

relations are given in Appendix [Cl Hence, in general, a given full-force mode ^f"n± carries contributions from 

tensorial-harmonic functions with ^ — 3 < Z' < ^ + 3. This coupling arises, of course, from the decomposition of 
the tensor-harmonic contributions into spherical harmonics. 
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E. Conservative and dissipative pieces of the SF 

In analyzing the physical consequences of the SF — as we start to do in Sec. E — it is useful to consider separately 
its conservative and dissipative pieces. We therefore now define these two pieces, and obtain a separate mode-sum 
formula for each. To this end, we first introduce the notation F" = F^^, reminding us that the SF is derived from 
the physical, retarded metric perturbation ha/3 = h'^afi ■ then define the "advanced" SF through 

00 

^^a"dv = E K± (Kt) -AIL-B-], (2.26) 
1=0 

where the modes -Ffuiii are constructed precisely as prescribed in the previous subsection, only this time using the 
multipole modes of the advanced metric perturbation, h'^^ . The regularization parameters A" and i?" are the same 

as those given above for the retarded SF. Since h'^^^ and have the same local singular behavior near the particle 
pol im , the sum in Eq. ()2.26|) is guaranteed to converge. 

Following Hinderer and Flanagan [2^, we define the SF's conservative and dissipative components as the parts of 
the SF which are (correspondingly) symmetric and anti-symmetric under retoadv: 

F"{= Fj"t) = + i^^igg, (2.27) 

where 

-^cons = 2 (-^rct + -^adv) J ^diss = 2 ("^rct ~ -^adv) • (2.28) 

Substituting from Eqs. (|2.1ip (with F" — ^ F^^.) and (|2.26p . we obtain the mode-sum formulas 

00 00 

'P'rans = E ['^full(cons)± ^ ^±-^ ^ -S"] = E -^regCcons) ; (2.29) 
/=0 1=0 



^diss — E ^Ml(diss)± = E '^rag(diss)> (2.30) 

where 

'^fun(cons)± = 2 ['^fun±(^'Q°/3) + '^fun±(^Q/J^)] : 'ffull(diss)± = 2 [-^l^\\±(^ap) ~ ^Mlil^a/T)] ■ (2-31) 

Notice that the dissipative piece of the SF requires no regularization within the mode-sum scheme. 

In Eqs. (|2.29p and (|2.30p the splitting of the SF into its conservative and dissipative pieces is performed mode-by- 
mode. This is useful in practice, because the Z-modes of the two pieces exhibit a rather different large-Z behaviour: 
While ^fun(cons)± normally admit an asymptotic power series in 1/Z [starting at 0{l)], the modes of -Ff^l^jj^^^^ die off 
at large I faster than any power of I. We will come back to this issue in Sec. IIII Dl 

The extraction of the conservative and dissipative pieces using Eqs. (|2.29p and (|2.30p entails a calculation of both 
retarded and advanced metric perturbations. This would normally double the computation time, as it requires one to 
solve the perturbation equations twice, changing the boundary conditions in order to obtain h'^p . Fortunately, in the 
case of a Schwarzschild background we can avoid this extra computational burden using a simple trick. For a given 
eccentric geodesic, we think of the SF as a function of r along the orbit. Without loss of generality we take t — 
to correspond to a certain periapsis passage [i.e., rp(T = 0) = r„iin]. Then we have the following symmetry relation, 
immediately following from Eq. (2.80) of [25]: 

KdAr) - Ho^)FZt{-r) (2.32) 

(no summation over a), where e((j) = (—1,1,1,-1) in Schwarzschild coordinates. This relation can be used to 
re-express Eqs. (|2.28p in terms of the retarded SF alone, in the form 

Fcl^Jr) = i [FZdr) + e^o.)FZt{-r)] , ^^^.(r) = ^ [F.^^ir) - e(„)i^."et(-r)] . (2.33) 

Similar relations are applicable to the I modes ^tun(cons)± ^^'^ '^fuii(cons)± ^® we^h These relations allow us to extract 
the conservative and dissipative pieces of the SF in practice without resorting to a calculation of the advanced 
perturbation: All that is required is knowledge of the physical (retarded) SF along the orbit. 63j 
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III. NUMERICAL METHOD 

In this section we describe the numerical method used to solve the field equations (|2.20p and to construct the 
local SF along the eccentric orbit. Our numerical scheme is basically similar to that of Paper I — we still use finite 
differencing on a characteristic grid in 1+1 dimensions — but we have modified our code in several important aspects. 
Most importantly, we abandon the use of a Richardson extrapolation over the grid size: this technique relies crucially 
on the uniformity of the local discretization error along the orbit, which can no longer be guaranteed in any simple 
way when dealing with eccentric orbits. To accelerate the numerical convergence we have instead upgraded our finite- 
difference scheme from second-order to fourth-order. This introduces a significant amount of additional complexity, 
especially in the treatment of grid cells traversed by the particle. Our method inherits from Lousto fS] (fourth-order 
scheme for l-flD evolution with a particle source) and Haas [12] (implementation for a scalar field), but we deviate 
from their methods in several aspects. 

A. Numerical domain and initial data 

Our integration domain is discreterized using a two-dimensional uniform mesh based on the double-null coordinates 
V and M, as depicted in Fig. [2l The numerical evolution starts with characteristic initial data 

h'-'^u = uo,v) = U'\u,v = vq) ^0 foraUi, (3.1) 

where the vertex (vq, uq) corresponds to the particle's location at < = [so vq ~ —uq ~ r^p{t — 0) where = ?'*(''p)]. 
The early stage of the evolution will be dominated by spurious radiation resulting from the imperfection of the initial 
data. However, as demonstrated in (2Q] (also in Paper I), these spurious waves damp down rapidly, and the error 
related to this behavior becomes negligible at late time. Our numerical algorithm monitors the residual error from 
spurious initial waves by comparing the SF values recorded at regular intervals along the orbit. We then make sure 
to evolve long enough for this error to drop below a set threshold. For a fractional error threshold of 10"'* in the final 
SF we find that the error from the spurious radiation can be safely ignored after ~ 2-3 x Tr of evolution (depending 
primarily on the value of p; larger p requires a longer evolution) . 




FIG. 2: Numerical domain: a staggered l-|-l-dimensional mesh in null coordinates v,u. r, is the standard Schwarzschild 
'tortoise' radial coordinate. The dotted line represents the trajectory of a typical eccentric orbit. In actual implementation the 
mesh is, of course, much finer than it is depicted here. 

Note that, in our setup, the numerical domain has no causal boundaries. Therefore, no boundary conditions need 
be imposed. 

B. Finite-difference scheme 

To derive our finite-difference equations, let us focus on a grid cell of dimensions Av x Au = h x h — say, the one 
in Fig. [3] with center C and vertices 1, 2, 3 and 4. We shall assume that the numerical values of /i*^*-' at points 2-15 
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are already known, and we wish to derive the value of h^'^^ at point 1. To this end we consider the integral of the field 
equations (|2.20p over the cell with center C. The ww-derivative term on the left-hand side is integrated in exact form 
to give 

/ U'Zdudv = h'^^^ -h^i^ -h'i^ +h'i\ (3.2) 

Jccll 

where hn^ denotes the value of ft-^'^ at the grid point labeled 'rt' in Fig. [S] The integral of the source term on the 
right-hand side of the field equations is expressed as 

cW ^ f <^('h(r-r ]dvdv - / 2/t*/ c?i/p-^5';^^(a;p(i)), (orbit crosses cell), 

"Li "" ^ rp)dz.di; - 1^ . (orbit outside cell), ^^'^^ 

where ti and tf are the values of t at which, correspondingly, the particle enters and leaves the cell in question. 
Since the integrand on the right-hand side depends only on the known trajectory of the particle (obtained in advance 
by solving the geodesic equation numerically), the integral can be evaluated in exact form. Our code implements a 
5-point closed Newton-Cotes formula ("Boole's rule") to evaluate this integral at each grid cell crossed by the particle. 




FIG. 3: Grid points involved in constructing our finite-difference scheme. The point C is the center of the cell over which we 
integrate the field equations, as described in the text. The dimensions of each grid cell are Av x Au = h x h. 



The rest of the terms appearing in the field equations (|2.20p [cf. Eqs. (|Bip - (jB10p ] can each be expressed schematically 
as either H = /(r)/i^''^'™, H ^ or H r, , where /(r) is some known function of r, j = 1, . . . , 10, and the indices {j)lm are 
suppressed for brevity. To complete the formulation of the finite-difference scheme we need to obtain finite-difference 
expressions for the integrals 

/= / Hdudv, ly = Hydudv, and Ir, = / Hr,dudv. (3-4) 

Jccll Jccll ' ' ./cell 

We are aiming here to achieve a quartic [0{h^)] global numerical error in the fields ft.'^'^'™. Hence, the local finite- 
diff'erence error for each of the integrals /, /„ and Ir, should not exceed 0{h^). To formulate the necessary finite- 
difference relations we will need to consider separately the following 3 cases (referring again to the grid cell with center 
C shown in Fig. [3]): (1) The orbit does not cross the triangular region shown in the figure ("vacuum cell"); (2) The 
orbit crosses either the segment 2-11 or the segment 3-15 ("near-orbit cell"); (3) The orbit crosses either the segment 
1-2 or the segment 1-3 ("orbit cell"). 



1. Vacuum cell 

Consider the formal two- variable Taylor expansion of a typical term H about the center of the cell in question — point 
C in Fig. [3l with coordinates v = Vc and u = Uc- We have 

N 

H{u,v)= ^,{u-u,r{v-v,)'> + 0{h''+'), (3.5) 

^ — ' a!o! 

a+b=a 

where a, b and A'' are non- negative integers (the latter to be specified below), and Cat are constant coefficients. Since 
the desired error in / (the integral of H over the 2D cell) is 0{h^), we are allowed an error of 0(/i^) in H and hence 
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take iV = 3 in Eq. p.5p . This leaves us with 10 expansion coefBcients Cab, which we can solve for, using Eq. (I3.5p . in 
terms of the values of H at the 10 points n = 1-10 indicated in Fig. [31 Substituting the values of these coefficients 
back in Eq. p.5|) and integrating over the grid cell, we obtain 

/ = — [2Hi + 10(^2 + H3 + Hi) - 4(i/5 + He) + {H7 - Hs - Hr> + Hio)] + 0{h'') [vacuum], (3.6) 

where iJ„ denotes the value of H at grid point n. To evaluate /„ and Ir, at local error 0{h^), we need instead truncate 
the Taylor series p.Sp at = 4, now leaving us with 15 coefficients Cab- These are determined by solving Eq. p.Sp 
given the values 7J„ at the 15 points n = 1-15. Taking dy and dr, in Eq. (j3.5p and integrating over the grid cell gives 

Iv - ^[9iH,-H2) + 19{Hs-H4)-5{He-Hg) + {H,o-Hu)] + 0{h^) [vacuum], (3.7) 
Ir, = ^[28iH3-H2)-5iHe-H5-Hg + Hs) + {Hw~Hr~Hu + Hi2)]+Oih'') [vacuum], (3.8) 
which are the desired finite-difference expressions for ly and Ir, ■ 



2. Near-orbit cell 




FIG. 4: Same as in Fig. [3l but now point C is located near the particle's worldline, represented by the dashed line. The 
finite-difference scheme for this case is described in the text under "near-orbit cell" . 

The above derivation assumes that the fields /i^*)'™ (and hence H) are sufRciently smooth across the region shown 
in Fig. [3l This is no longer the case if the worldline of the particle traverses this region, since H is generally non- 
differentiable across the worldline, and the above Taylor-series-based method would fail. Let us first consider the 
simpler case, where the particle crosses either the segment 2-11 or the segment 3-15 (i.e., it does not traverse the 
integration cell itself), as demonstrated in Fig.^l In this case the function H{u,v) can still be expressed as a formal 
Taylor series, 

^ r± 

H{u,v)^ E ^(^""c)"(^^-«c)' + 0(/i('^+'^ (3.9) 
^-^ a\b\ 

a+b=a 

but we now have two different sets of expansion coefficients, c^^ and c^f^, depending on whether r{u,v) > rp{t) 
or r{u,v) < rp{t), respectively. As in the vacuum cell case, the value of H at the grid points 1-15 provides 15 
independent equations for the unknown coefficients c^^^, which, however, are now 30 in number. The necessary 
additional 15 relations between the various c^{,'s are obtained by utilizing explicit junction conditions for H and its 
derivatives across the particle's orbit, as we now describe. 

In Appendix |E] we derive explicit jump conditions for the perturbation /i^*)'™ and its derivatives at a generic point 
xq along the (known) geodesic worldline. Specifically (and using the notation of Appendix |E]) , we calculate the two 

jumps [h^u^"^]o and [h^v''"^]o, as well as the three jumps [h^uu^]o, [^!'™™]o and [/i[™'"]o, the 4 jumps in the various third 
derivatives, and the 5 jumps in the various fourth derivatives. Together with the continuity condition [/i*^*-"™]o = 0, 
we hence obtain 15 jump conditions in total. Now referring back to our near-orbit cell scenario and to Fig. SI we take 
xq to be the intersection of the worldline with the past light cone of point 1, as demonstrated in the figure. The 15 
jump conditions for and its derivatives at xq readily translate into 15 jump conditions for H and its derivatives 
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at that point. Imposing these conditions in Eq. (j3.9l) yields the required additional 15 independent linear equations 
for the coefficients c^^^. Our algorithm solves the total of 30 equations for c^^^ numerically, given the numerical values 
of Hn at points 1-15. 

can be integrated over the cell of centre C, giving 



Once the coefficients c^^ have been calculated, Eq 



I = h' 

Iv = 



-00 ■ 



24 



'^02) 



^oi + 24 ("^^i + 



0{h^) [near-orbit], 
^0{h^) [near-orbit]. 



-10 



^ (^21 - 3c: 



30 



3c; 



03 



-fO(/i^) [near-orbit]. 



(3.10) 
(3.11) 
(3.12) 



where the values c^j^ apply if C is located at r > rp(i) and the values c^^ apply if C lies in the region r < rp{t). 



3. Orbit cell 



The procedure described above for calculating the Taylor coefficients c^^ is applicable even when the particle's 
worldline crosses the considered grid cell. However, the integration of H and its derivatives over the grid cell then 
becomes slightly more involved, because the cell is divided by the trajectory into two parts, in each of which the 
coefficients take different values. Following Lousto [42], we consider separately the four cases illustrated in Fig. [51 




FIG. 5: Illustration of the four cases considered in formulating the finite-difference equation for grid cells traversed by the 
particle's worldline (here represented by dashed lines). The center of the cell is at {u, v) — {uc,Vc), and in each case we indicate 
the u,v coordinates of the two points where the particle enters/leaves the cell. 

"Case UU" (top left in Fig. [5]): the particle enters the cell crossing the u=const segment 2-4, and leaves it crossing 
the w=const segment 1-3. The worldline splits the cell into two bits, "left" and "right", respectively labeled L and 
R in the figure. Denoting the corresponding contributions to / by and I^, we have I ^ + , where, using the 
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expansion p.9|) . 



rh/2 


fh/2 


/ dv 


/ duH = 


J-h/2 


/gp(S) 


fh/2 


/•"p(«) 


/ dv 


/ dui/ = 


l-h/2 


'-h/2 



a+b=0 
3 



(a + l)!fe! 



(a + l)!6! 



a+6=0 



/i/2 
-/i/2 

/!./2 



dvv^ 



(3.13) 



h/2 



{u^{v)Y+' -{-hl2f+' +0{h^), (3.14) 



Here u = u — Uc and w = w — Uc, and 'Sp(u) represents the value of u on the trajectory, viewed as a function of v. 
Similar expressions are easily obtained for /„ and Ir, ■ 

"Case VU" (top right in Fig.[S]): the particle enters the cell crossing the u=const segment 3-4, and leaves it crossing 
the ?;=const segment 1-3. We denote the entry w-value by Vi and the exit u-value by u/, and further denote Vi = Vi — Vc 
and Uf — Uf — Uc- We obtain, in this case, 

h/2 i-h/2 

dv I duH 

-h/2 
3 



a+b=0 



(a + l)!(6+l)! 







[ ^ du 1 


dvH 


J-h/2 J 


-h/2 







+(a + l) / duu'' {v^{u)f^^ - {-hj-i) 

-h/2 L 



b+1 



0(/i«), 



du 



h/2 



dvH^ 



-h/2 Jv^iu) atbio^-^^^^^- -^-'^Z 



duu°- 



{hl2f+'-{v^{u)r' +Oih'), 



b+1 



(3.15) 
(3.16) 



where Vp{u) is the value of v on the trajectory, expressed as a function of u. Once again, similar expressions can be 
obtained for ly and /, _^ . 

"Case VV" (bottom left in Fig. [S]): the particle enters the cell crossing the u=const segment 3-4, and leaves it 
crossing the u=const segment 1-2. In this case we obtain 



h/2 fVpiu) 

du 

-h/2 J-h/2 
h/2 ph/2 



a+b=0 



a!(6+l)! J_ 



h/2 



du M° 



du 



dvH^ ^ 



-af) 



h/2 Jv^iu) afblo "'(^^ •^-'»/ 



h/2 
h/2 



{vp{u)f+' - (-V2) 



dun"- 



\h/2r' - {v,{u)r' 



b+l 



(3.17) 
(3.18) 



and similar expressions for /„ and 7^, . 

"Case UV" (bottom right in Fig. [5]): the particle enters the cell crossing the w^const segment 2-4, and leaves it 
crossing the u=const segment 1-2. We denote the entry u- value by Ui and the exit u- value by v/, with Ui — Ui — Uc 
and Vf = Vf — Vc- In this final case we have 



dv 



h/2 



du 



-h/2 Jup{v) 
h/2 j.h/2 

dv du H 

f J-h/2 



a + l 



a+b=0 
f 



[a + iy.bi 

p(C) 

duH 



dvv" 



{h/2r+' - {up{v)) 



a+l 



dv 

-h/2 J-h/2 



y ^ 



a+b=0 



(a + l)!(5+l)! 



{h/2)''+^ [1 + (-1)"] {h/2f+^ - {vff^^ 



-{h+l) f dvv" {up{v)T^' - {-h/2y 

-h/2 L 



(3.19) 



(3.20) 



with similar expressions for /.„ and Ir 
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4- Predictor-corrector method 
In summary, recalling Eq. (j3.2p and with reference to Fig. [3l our basic finite-difference formula takes the form 



dudv 



coll 



(3.21) 



where iS'*-* is given in Eq. p.Sp and the integral (over the cell with center C in the figure) is evaluated with local error 
0{h^) as discussed above. A complication arises, since our finite-difference expressions for the integral in Eq. (I3.2ip 
involve the value of the perturbation at point 1 itself [as in, e.g., Eq. p.6p ]. which is the very unknown value we wish 
to compute. 

We overcome this difficulty using a type of predictor-corrector algorithm, whereby we first approximate the value 
of the field at the point in question (our point 1) using extrapolation, and then apply our finite-difference formula 
(I3.21[) iteratively, until the required accuracy is achieved. Specifically, we use the values of the perturbation H at the 
four points 3, 6, 10, and 15 (see Fig. [S]) to extrapolate the value Hi with an error of 0{h'^). If the particle's worldline 
happens to cross the null segment 1-15 we instead use the points 2, 5, 7, and 11 for this extrapolation. We then use 
the value of Hi thus obtained as input in Eq. p.2ip . Since the extrapolated terms enter Eq. p.21|) multiplied by at 
least one power of h [see, e.g., Eqs. p.6p and p.7p ]. the resulting value of Hi would have a local error of 0{h^). We 
then apply Eq. (j3.21l) once more, with the new value of Hi as input. The output value of Hi following this second 
iteration should now have a local error of 0{h^) as desired. 

The above finite-difference scheme is designed to yield a local error of 0{h^) in hi^ at each grid point. Since the 



overall number of grid points contributing to the accumulated error scales as 
quartic [i.e., 0{h^)] numerical convergence. 



h , we expect our scheme to show a 



C. Monopole and dipole modes 

In the monopole case {I = 0) the system of 10 field equations (I2.20p reduces to 4 equations only; it reduces to 6 
equations for each of the two dipole modes {l,m) ~ (1, ±1) and to a single equation for the dipole mode {l,m) = (1,0). 
One may attempt to apply the above numerical evolution scheme to these non-radiative modes as well. However, 
numerical experimentation suggests to us that the monopole and dipole cannot be evolved stably using this scheme. 
A naive application of the evolution scheme yields exponentially growing solutions, and, since our scheme gives us no 
handle on the boundary conditions, the occurrence of these unphysical solutions is difficult to control. 

Instead, we deal with the two modes / = 0, 1 using the standard frequency-domain method, just as in Paper I. The 
physical Lorenz-gauge monolpole and dipole are constructed from a basis of homogeneous frequency-mode solutions 
of the underlying ordinary differential equations. These homogeneous solutions are obtained numerically. However, 
unlike in the circular-orbit case dealt with in Paper I, here we face the "Gibbs phenomenon" , since for an eccentric 
orbit the perturbation is a non-differentiable function of coordinate time t at the particle's location. We have discussed 
this problem in depth in Ref. [3l| . and proposed an elegant solution, whereby the correct physical perturbation is 
constructed as a sum of "fake" frequency- mode solutions whose Fourier sum converges exponentially even at the 
particle's position. Here we apply this method in order to obtain the physical Lorenz-gauge monopole and dipole 
modes. A full description of this construction will be given in a forthcoming paper [43j . 



D. Implementation of the mode-sum scheme 

Once we obtain the numerical values of the functions /i'*^ and their derivatives along the orbit (over a complete 
radial period), we can construct the full- force modes -Ff"n-|- at any point along the orbit through the procedure 
described in Sec. Ill Dl The mode-sum formula (|2.1ip then gives the physical SF at that point. The application of the 
mode-sum formula involves summation over contributions from an infinite number of modes, from I — to I — oo. 
In reality, of course, we are only able to compute a small number of low multipole modes — not least because the 
numerical calculation becomes increasingly more demanding with larger Since the mode-sum scheme converges 
rather slowly (as ~ 1/0 j a calculation of the SF with even a modest accuracy requires that we take a proper account 
of the contribution from the truncated tail of the mode sum. 

Let us denote by I the highest spherical-harmonic Z-mode calculated numerically (note this would entail calculating 
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all tensor-harmonic modes from Z = to Z = / + 3). Recalling the notation of Eq. (|2.1ip . we express the SF as 

/ oo 

^" = E K + E ^ ^ili + (3-22) 
(=0 ;=/+! 

where F^^^ is the part calculated numerically, and F^j is the truncated tail we need to estimate. Here it is beneficial 
to consider separately the conservative and dissipative pieces of the SF. We remind that these can be constructed 
individually using the mode-sum formulas (j2.29p and (|2.30p via the procedure described in Sec. Ill El 

Consider the conservative piece first. The regularized force modes in this case admit the large-Z expansion 

^"eg(cons) = + D%L-^ + . . . , (3.23) 

where D'^2n /-independent coefficients. An approximation for these coefficients can be obtained by fitting a large-Z 
subset of numerical data to Eq. p.23p . In practice we take I = 15 and find the two coefficients ^-nd £'"4 using 
the numerically-derived modes 10 < ^ < 15. The large-Z tail piece of ^'regCcons) ^'^^^ approximated as 

00 

Kons,i>i ^ E ^"eg(cons) ^ Dt^T i{T + 3/2) + D%T:,{1 + 3/2)/3!, (3.24) 
1=1+1 

where r„(a;) is the polygamma function of order n, defined in terms of the derivatives of the standard gamma function 
as r„(x) = (i"+^[logr(x)]/(ix"+^. Since the leading term omitted in Eq. p. 231) is of 0(L~^), we expect the absolute 
error in our estimation of F^^^^^^j to be of 0{l~^), or 0{l~'^) fractionally. With / — 15 this amounts to a ^ 10~^ 
fractional error, which we can afford to tolerate in this work. 

Now consider the dissipative piece. We have that the magnitude of ^'i^gjdiss) ^^^^^ faster than any power of 1/Z 
at large I. For the range of orbital parameters explored in this work we find that the numerical value of F"^^^^^^^-^ 
drops below the round-off error at I ^ 7-12, and so the large-Z tail, estimated as F" '-^ F^K.- s, can be safely 
neglected taking I = 12. In practice, to avoid adding up spurious round-off contributions, we truncate the mode-sum 
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series at I = min{/, 15}, where I is the first value of / above I = 7 for which we find 

Our procedure for estimating the numerical error is as follows. First, we estimate the discretization error in each 
of the computed /-mode contributions (as a time series along the orbit) by repeating our numerical evolution at a 
coarser resolution and using the difference between the high-and low-resolution data sets as a crude error estimator. 
[For example, to obtain the data in Tables Ill IIVI below we applied a cell size of Au x Av — (0.02M)^ for our highest 
resolution runs, then Au x Aw = (0.04M)^ for error estimation.] The total error in the numerically-computed part 
of the SF is then taken (conservatively) as the sum of the errors (in absolute value) from the various modes. In the 
case of the conservative piece we add to this the estimated standard fitting error for the large-/ tail. For the fitting 
itself we use the high- resolution /-mode data points, weighted by their estimated discretization errors. This procedure 
yields a conservative estimate for the numerical error in each of i^c"ns(x) F^^^^{x)- It is this error that we quote 
in the next section when presenting our results. 



E. Code validation and performance 



Using a few test runs with a range of sample parameters, we tested our code (i) for 4th-order numerical convergence; 
(ii) against the SF results of Paper I in the circular-orbit case (e = 0); (iii) by extracting the flux of energy and angular 
momentum carried by the gravitational waves and comparing with results in the literature (see Sec. IIV A]) : and (iv) 
by verifying that the dissipative component of the computed SF precisely balances the above flux (Sec. IIV A[) . Our 
code seems to perform well, at a standard fractional accuracy of < 10^^ in the final SF, across the parameter range 
< e < 0.5 and p < 20M. For larger eccentricities and/or larger values of p the long evolution time required begins 
to play a prohibitive role. Our code is still fully functional at (say) {p, e) — (50, 0.9), although in this case it becomes 
computationally impractical to achieve our standard lO"** accuracy running on a standard (single) desktop computer. 

In developing and testing our code we used a desktop workstation with a 2.5GHz 64-bit processor and 8Gb of RAM. 
A typical computation of the SF over a complete radial period, with given parameters in the above "workable" range 
(and at the above accuracy standard), demands 4-8 days of CPU time on this machine. 
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IV. SAMPLE RESULTS 



Figure [6] displays SF results for the sample parameters (p, e) = (10,0.2), (10,0.5) and (10,0.7). We plot the 
temporal and radial components of the (total) SF along the geodesic orbit as functions of coordinate time t. The 
third non-trivial component, F"^, can be obtained using the orthogonality condition F°'Ua = 0. In these plots t = 
corresponds to a periapsis passage (where r = Tmin), so that, in accordance with the discussion at the end of Sec. Ill El 
the dissipative/conservative pieces of are described by the symmetric/anti-symmetric parts of the i<"*-graph, and 
conversely for F^ . 
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FIG. 6: The gravitational SF for a sample of eccentric orbits. We plot the temporal and radial components of the SF 
as functions of t for the sample parameters (p, e) = (10,0.2), (10,0.5) and (10,0.7). The solid and dashed lines show 
and , respectively. The (f component of the SF can be trivially obtained using the orthogonality condition F"ua = 0. 
In all graphs t = corresponds to a periapsis passage. The radial periods for the e = 0.2,0.5,0.7 plots are, respectively, 
Tr — 328M, 434M, 693M. We have cut out from these plots the early part of the numerical solution, where non-physical 
initial spurious waves dominate; the plots display only the later, stationary part of the solutions, after the spurious waves have 
dissipated away. Note the small retardation manifest in the amplitude of the total SF (with reference to the orbital phase). A 
similar retardation is observed in the scalar and EM cases [l^ . v2^ . 

In Tables HHTVl we present numerical SF results for the sets of orbital parameters (p,e) = (7,0.2), (7,0.4), (10,0.3) 
and (15,0.3), for a few sample values of the radial phase x- These numerical values can be used as a reference 
for testing future calculations of the SF. In the tables we display the conservative and dissipative pieces of the SF 
separately. We only display values for F^ and F^ — once again, the azimuthal component can be obtained trivially 
using F°'Ua — 0. Also, we only give values for the "outbound" half of the orbital period (0 < x < ti"; where Tp > 0); 
the values for the "inbound" half (tt < x < 27r) can be obtained immediately based on the symmetry relations ()2.33|) . 
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TABLE I: The gravitational SF at selected points along an eccentric geodesic with (p, e) = (7,0.2). rj = fi/M. In the first 
column X is the radial phase along the orbit [cf. Eq. (|2.5|l ]. with x = 0, tt corresponding to periapsis and apoaspsis, respectively. 
Subsequent columns display, in order, the conservative and dissipative pieces of F* and the conservative and dissipative pieces 
of . Values in brackets are estimates of the uncertainty (due to numerical error) in the last displayed decimal place. For 
example, 5.3844(3) x 10-* stands for (5.3844 ± 0.0003) x 10-*. The values for F"^ can be obtained from F^-Uc = 0. Values for 
the SF along the "inbound" half of the radial period (tt < x ^ 27r) can be deduced based on the symmetry relations (|2.33|) . 
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TABLE II: The gravitational SF at selected points along an eccentric geodesic with (p, e) = (7, 0.4). The structure of the table 
is the same as that of Table. H] 

A. Dissipation of energy and angular momentum 

Given the local SF, we can calculate the (orbit-averaged) rate at which orbital energy and angular momentum are 
dissipated. This information is contained in the t and tp components of the local SF. From Eq. (|2.10p one readily 
obtains 

£{X) - -WixT'Ftix). t{x) = WixT'F^ix), (4.1) 
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TABLE III: The gravitational SF at selected points along an eccentric geodesic with (p, e) = (10,0.3). The structure of the 
table is the same as that of Table. IT] 



18 



X 


77-2 

'/ ' cons 




77-2 F'' 

'I cons 










-1.040267(9) X 10"* 


1.139648(3) X 10"2 





7r/8 


3.03529(3) X 10"'' 


-8.23882(8) x 10"^ 


1.102460(1) X 10"^ 


1.03215(6) X 10"" 


7r/4 


5.07593(2) X 10"* 


-3.9157(1) X 10-5 


1.000022(2) X 10"^ 


1.43970(1) X 10"" 


37r/8 


5.68296(2) x 10"" 


-8.3320(6) X 10"^ 


8.56001(2) X 10"^ 


1.19734(3) X 10"" 


7r/2 


5.06147(1) X 10"'' 


1.97928(7) X 10"'' 


6.99810(1) X 10"^ 


7.2183(3) X 10"5 


57r/8 


3.765543(7) x 10""' 


2.2286(2) X 10"® 


5.578997(8) x 10"=* 


3.4786(1) X 10"5 


37r/4 


2.340696(1) X 10"" 


8.4719(9) X 10"^ 


4.486884(2) x 10"^ 


1.4418(1) X 10"5 


77r/8 


1.0842785(8) x 10"'' 


1.0087(2) X 10"^ 


3.814949(2) x 10"^ 


4.9887(8) X 10"** 


TV 





-9.092(5) X 10"** 


3.590240(2) x 10"^ 






TABLE IV: The gravitational SF at selected points along an eccentric geodesic with (p, e) = (15,0.3). The structure of the 
table is the same as that of Table. HI 



where in this section (unlike elsewhere in this work) an overdot denotes d/dt. We shall assume here, in effect, that 
fi/M is sufficiently small that the back- reaction effect on the orbit over a period of a few Tr can be neglected. At this 
"adiabatic" limit the functions £{x) ^'^id C{x) B^re periodic with period Tr, and their time- average is hence given by 



1 



SdX: 



1 

Tr Jo 



Cdx. 



(4.2) 



These averages determine the "secular" dissipative drift in the values of £ and C. Note that {£) and (£) depend only 
on the dissipative components F^^^^ and F^^^^ (respectively); it is immediately evident from the symmetry relations 
(|2.33p that the contributions from F^:^^^ and F^f^^^^ vanish upon orbital-averaging. We also note that the dissipative 
radial component F^^^^ (let alone i^cons) effect on the values of {£) and (£). The secular drifts {£) and (£) must 

be balanced by the flux of energy and azimuthal angular momentum in the gravitational waves radiated to infinity 
and down the event horizon. Denoting the respective energy fluxes by (£')oo/eh ^-nd angular momentum fluxes by 
(L)oo/EHj we have the balance equations 



-^,{£) = {E)^ + {E)BH = {E)tote.h 
-H{C) = (L)oo + (i)BH = (i)total- 



(4.3) 
(4.4) 



Two validation tests for our code now suggest themselves. First, we may attempt to extract the asymptotic fluxes 
directly from our numerically-calculated Lorenz-gauge metric perturbation, and compare with results in the literature. 
Second, using Eqs. (|4.ip we can derive the local quantities {£) and (£) from our SF results, and check whether the 
balance equations (|4.3p and (|4.4p are indeed satisfied. 

For both above tests we need a time-domain formulation of the asymptotic fluxes in Schwarzschild spacetime. 
Such a formulation was presented by Martel [3] and Poisson [4^ 64 1, and we shall adopt it here. Martel-Poisson's 
construction is based on the Regge- Wheeler and Zerilli-Moncrief perturbations functions ^'({^ 
related to our Lorenz-gauge variables through 



and ^zM' which are 
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Vl/"" — 
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2(; + 2)! _ 
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^(9) + Lhim _ 

r r 
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2r ' * 2Ar 



(4.5) 
(4.6) 



with \ = {I + 2){l - 1). In terms of and *z'Si, the fluxes at infinity are given by [4J, |4£ 
(^)- = ^Ef-^(4|*Rw(")P + I^ZM(")P 



Im 



{L)c 



-E 



647r ^ (/ - 2)! 

Im 



(4.7) 
(4.8) 



where an asterisk denotes complex conjugation 
a period Tr would suffice 



and the functions ^'[^^ ^'^d evaluated at the 



indicates a suitable time-average (in our case an average over 

wave zone" , w — 00 with fixed u. 
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This paper (TD) 


Martel (TD) 


Fujita et al. (FD) 


p = 7.50478, 


e = 0.188917 










-2 


3.169(1) 


3.1770 


3.16899989184 


{L)oo X 10^7?- 




5.96760(8) 


5.9329 


5.96755215608 


p = 8.75455, 


e = 0.764124 








{E)^ X loS 


-2 


2.124(3) 


2.1484 


2.12360313326 


(i)oo X 10^77- 




2.7774(6) 


2.7932 


2.77735938996 



TABLE V: The flux of energy and angular momentum in the gravitational waves radiated to infinity: comparison with results 
in the literature. The second column shows the values of the radiative fluxes (-B)oo and {L)oo, evaluated from our numerical 
results using Eq. (|4.7p and (|4.8|l for two sample values of p, e. Values in parentheses estimate the uncertainty in the last 
displayed figure. The subsequent columns display, for comparison, the corresponding values obtained by Martel !44| and Fujita 
et al. 47]. TD/FD indicate time/frequency-domain methods. Fujita et al. claim all their displayed figures are significant. 

The horizon fluxes (-E)bh and {L)bh are given by expressions similar to (|4.7p and (|4.8p . merely replacing u ^ v and 
evaluating ^'[^^ ^^"^ ^zm near the horizon, i.e., at w — > 00 with fixed v. 

To calculate {E)oo and (i)oo, we start by recording the numerical values of the perturbation functions h^^\u) at 
V = lOOOOAf (approximating null infinity), over a complete radial period. The desired fluxes at infinity are then 
calculated using Eqs. (I4.5p - (l4.8p . making sure that the error from truncating the sum over / is properly controlled 
(this is not difficult, as the mode sum converges exponentially; in none of the cases considered here we found it 
necessary to include modes beyond I = 12). In a similar manner we also construct the horizon fluxes (£')eh and 
(L)eh, starting by extracting the numerical values h!'^\v) at very large u (approximating the horizon; in practice we 
used u = lOOOOM), and then using Eq. (|4.5p and (|4.6p with the horizon- version of Eqs. (|4.7I) and (|4.8I) . We estimate 
the error in our flux values by comparing results obtained at two different grid resolutions {h = 0.1 against h = 0.2). 

The eccentric-orbit fluxes {E)oo and (L)oo were computed independently in the past by several authors, including 
Tanaka et al. (Sand Cutler et al. [11] [using frequency-domain (FD) analyses based on Teukolsky's formalism], and 
later by Martel 4J] [using a time-domain (TD) analysis based on the Regge- Wheeler- Zerilli formalism]. Very recently, 
Fujita et al. [47] developed a highly accurate FD algorithm for flux calculations. In Table fVl we compare our flux 
data with the TD data of Martel and the FD data of Fujita et al. [i^. We look at two strong-field orbits, one with 
moderate eccentricity (e ~ 0.19) and the other with a rather high eccentricity (e ~ 0.76). All the results shown agree 
with ours to within 1%. The FD results agree with ours to within at least 0.01%, and in all cases the FD results fall 
well within our estimated error bars. Martel's TD results were presented without error bars, but they are likely less 
accurate than the FD ones. 

In Tables IVll and IVTIl we carry out the second test mentioned above, i.e., we check whether our numerical SF data 
and flux data satisfy the balance equations (|4.3p and (|4.4p . The tables compare the time-averaged rates of loss of 
orbital energy and angular momentum, /J.(f ) and /i(/C), with the corresponding total (horizon-|-infinity) fluxes (£^)totai 
and (L) total- We find that the two are consistent with each other to within the numerical accuracy. This constitutes 
a significant, highly non-trivial validation test for our code. 

Tables IVTl and IVTIl also display the partial contributions to the fluxes (i?) total and (L) total coming form black hole 
absorption, i.e., (ij)BH and (i)BH- Note that our numerical accuracy is sufficient to confidently resolve these (relatively 
minute) horizon fluxes in the examples considered. The consistency of the local dissipative SF with the asymptotic 
fluxes, at the accuracy level maintained here, is evident only when black hole absorption is correctly accounted for. 
The data in Tables IVII and I VIII represent a first numerical test of the Martel-Poisson TD horizon- flux formula. 

B. Zoom-whirl orbits 

An interesting family of eccentric geodesies, so called "zoom- whirl" orbits [1^, has 6 — 2e = e<Cl. These 
geodesies correspond to points in the p, e plane lying very close to the separatrix (see Fig.[T]), possessing energy-squared 
5^ only slightly smaller than the maximum of the effective potential i?(r, C^). A particle on a zoom- whirl orbit spends 
most of the radial period "whirling" around the central hole in a nearly circular orbit near periapsis, before "zooming 
out" back to apoapsis distance. During the whirl episode the particle may complete many revolutions in — the 
(/?-phase accumulated over one radial period scales as (x ln(64e/e) [see Eq. (2.25) in [ssj]. 

In Fig.[7|we show the SF along a sample zoom-whirl orbit with parameters (p, e) = (7.4001, 0.7). All components of 
the SF are relatively very small near the apoapsis, as would be expected in virtue of the large distance from the central 
black hole. During the brief "zoom-in" and "zoom out" episodes the particle has a large radial velocity component. 
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TABLE VI: Testing our code using the global energy balance equation (|4.3p . The third column displays the (negative of the) 
average rate of loss of orbital energy over one radial period, —fj,{£}, for strong-field orbits with semi-latus rectum p = 7 and a 
range of eccentricities e. This quantity is calculated from local SF data using Eqs. I|4.1|) and (|4.2|) . For comparison, we give in 
the fourth column the corresponding values of the total energy fluxes {-E)totai in the gravitational waves radiated to infinity and 
down the event horizon. These fluxes are extracted from our numerically-derived Lorenz-gauge metric perturbation, evaluated 
at the corresponding asymptotic domains. Values in parentheses are estimates of the numerical error in the last displayed figure. 
The last column shows the relative difference between —fJ^{£) and (E) total, confirming that the balance equations are satisfied to 
within our working precision. For reference, the fifth column shows the relative contribution to {i5)totai from horizon absorption 
alone, denoted {E)bii- Manifestly, for the strong-field orbits considered here and with our working precision of < 10"^, black 
hole absorption cannot be neglected. 
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TABLE VIL Testing our code using the global angular- momentum balance equation (|4.4p . The structure of this table is similar 
to that of Table IVTl The table compares between the average rate of change of orbital angular momentum, {£), as inferred 
from local SF data, and the flux of angular momentum carried away by the gravitational waves, {i/)totai, as inferred from the 
asymptotic waveforms. We also indicate the relative contribution to {L)totai from angular momentum absorbed by the black 
hole, denoted (L)bh. 

and the SF changes rapidly. During the whirl phase the particle moves on a nearly-circular orbit, and we expect 
the SF to settle to a constant value (this expectation is indeed confirmed in studies of the scalar-field and EM SFs 
[l^,!!^). We find that while this is true for the t and (p components, the r-component of the SF shows an unexpected 
linear-in-t behavior during the whirl. Other zoom-whirl orbits we examined showed a similar behavior. 

In order to understand the above peculiarity we conducted several experiments. The following points summarize 
the information gained, (i) The linearly growing piece of is entirely dissipative; the conservative piece -Fc'ons shows 
no such linear growth (see Fig. [8]). (ii) The linearly growing piece of F^-^^ is entirely attributed to the time-dependent 
piece of the monopole (l = 0) contribution to the SF; other modes show no such behavior, (iii) When examining a 
sequence of zoom-whirl orbits approaching the separatrix (e — >■ 0) with fixed e — see Fig. |5] — we observe that the rate 
of linear growth decreases with e, although rather slowly (slower than oc e). It is possible that our numerical results 
are in fact consistent with F^-^^ — >• as e — ?> (as expected), but this is difficult to verify numerically, (iv) We have 
applied our numerical algorithm to study the zoom- whirl behavior of the scalar-field SF. As in we observed no 
linear behavior during the whirl. ^65*] 

The last point (iv) serves to reaffirm our trust in the numerical code. The combination of points (i) and (ii) (together 
with the fact that no linear growth is observed in the scalar and EM cases) implies that the linear behavior is purely 
gauge-related, and suggests that it should have no observational (gauge invariant) consequences. In particular, the 
culprit monopole piece of F^-^^^ is non-radiative and has no secular physical effect on the orbit. Finally, the observation 
made in (iii) suggests there is nothing wrong with our choice of gauge either. In fact, it may be argued that the observed 
linear-in-t behavior, with a weak dependence on e, is perfectly consistent with the theoretical expectation based on a 
local analysis of F^^^^ near the separatrix. We explain this in the following. 

Consider the behavior of F^-^^^ for a zoom- whirl orbit, e ^ 1. During the whirl the radial phase x changes very 
little, so, taking x = at periapsis as usual, we can assume x 1- In the following analysis we fix the eccentricity 
e(> 0) and consider the limit e by taking p Q + 2e (from above), inspecting the behavior of F^iss leading 
order in both x and e. A convenient starting point for this analysis is the orthogonality condition u^F'^ = 0, whose 
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FIG. 7: The gravitational SF for a zoom-whirl orbit with parameters (p, e) = (7.4001,0.7). In the upper panel we show the 
(contravariant) t, r, and components of the SF as functions of time t along the orbit. The lower panel shows the radial 
motion of the particle, for reference, t = is periapsis, and the radial period is ~ 789M. While _F* and F*^ quickly settle to a 
constant value during the quasi-circular whirl episode (as one would expect), the radial component exhibits a peculiar linear 
behavior. This behavior is analyzed in the text (cf. Fig. [5]). 




FIG. 8: Conservative and dissipative pieces of the radial SF for zoom-whirl orbits — left and right panels, respectively. Shown 
are results for four orbits with the same eccentricity (e = 0.7) but with increasing proximity to the separatrix p = 6 -I- 2e. 



dissipative part can be rearranged to give 

-Ipr ^ (gV/p)(gdiss-l^/:diss) ^ {S^h)-^ .^g. 
^ diss ' ^ ' ' 

Here we have used the dissipative part of Eqs. (|4.ip . and denoted VI = d(pp/dt. Recall /p = 1 — 2M/r-p and an 
overdot denotes d/dt. The factor (£^//p) is regular at x = e = and hence uninteresting. For the radial velocity, Eq. 
(|2.8p gives oc e^^^x, where throughout our present discussion a proportionately symbol implies the leading-order 
scaling with x ^^nd e. The function ^' is symmetric in x [recall Eq. (|2.33p ] and clearly 5'(x = 0) = 0; hence we write 
^ — X^^(e, e) + O(x^), where ^ is x-independent. Taking now t(x = 0) = 0, Eq. (12. 6p gives x e^/^f, and we thus 
obtain cx et and '3/ cx et'^ip{e, e). We conclude that 

F,^i,,cxtV'(e,6). (4.10) 

It is now essential to understand the scaling of •0 with e as e — >■ 0. Clearly, ■0 — at this limit, since for circular orbits 
we have fdiss = ^^-Cdiss- It may be argued [38l. |49| that ip scales like the (small) fraction of the radial period that the 
particle is spending in the zoom phase, which, in turn, is proportional to 27r/Aiy9 oc [ln(64e/e)]^^. If this argument is 
to be trusted, we obtain 

F^,,, « t X [ln(64e/e)]-\ (4.11) 
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which may explain the very weak e-dependence of the slope in Fig. [8l It is difficult in practice to test our numerical 
results more quantitatively against the scaling relation (j4.11l) (precisely because the e-dependence is so weak), but we 
cannot rule out the possibility that our results are in fact consistent with this scaling. 

It remains to understand why the linear mode does not exhibit itself so pronouncedly in the scalar and EM cases, 
and how the amplitude of this mode might depend on the choice of gauge in the gravitational case. This requires 
further analysis, which we do not attempt here. We remind that, from a practical point of view, the linear mode 
should not cause any real concern, as it cannot affect any gauge-invariant quantity derived from the SF. 



V. ISCO SHIFT 



Using our SF code we can start to explore the 0(/i) "post-geodesic" dynamics of the orbit, and quantify the physical 
effects resulting from the finiteness of /j,. As a first concrete application of the code, we calculated the 0(/x) shift in 
the location and frequency of the ISCO due to the conservative piece of the SF. We recently reported the results of 
this calculation in a Letter (ssj. Here (and in Appendix IF| we provide full details of this analysis. 

The radiative transition across the ISCO, from a slow quasi-circular inspiral to a rapid plunge, has been studied 
in the past and is now well understood. The transition occurs not at a well-defined radius but rather in a gradual 
manner, across a "transition regime" , whose width (in terms of the corresponding azimuthal frequency) is proportional 
to a low power of the mass ratio: Afitrans c>c (/x/M)^/^ [s^, [HI- However, apart from the dominant radiative effect 
which drives the inspiral (and the eventual transition to plunge), the gravitational SF also has a conservative effect, 
which shifts the location of the ISCO away from r = 6M [by an amount of 0{ii)]. Unlike the radiative transition, this 
conservative shift is precisely quantifiable. Moreover, the value of the azimuthal frequency at the shifted ISCO, riiscoi 
is essentially gauge invariant, and hence provides a useful handle on the strong-field conservative dynamics. Indeed, 
the value of ilisco (for mass ratios not necessarily extreme) has long been utilized in testing and calibrating various 
approximate treatments of the general- relativistic binary problem (see., e.g., d^Hl], and the very recent [2l|). Our 
SF code allows us, for the first time, to obtain a precise value for ilisco (modulo a controlled numerical error) at 0{fi). 

In what follows we first derive a formula for the ISCO frequency riisco including 0{n) conservative SF corrections, 
and then describe the numerical method used to obtain the necessary SF data, and the results. Many of the details 
are relegated to Appendix [Fj 



A. Formulation 



We first review the notion of ISCO in the unperturbed [geodesic, 0{fjP)] case. The radial geodesic equation [obtained 
by differentiating Eq. (j2.2|) with respect to r] reads 



dr2 



-^cff(r-p,£^); Fcsir.C^ 



2 ■ 



(5.1) 



We consider a slightly eccentric orbit representing an e-perturbation of a circular orbit with radius Tq. We write 

rp(T) =ro + eri(T)+0(e2), (5.2) 
where e ^ 1 and ri(T) is e- independent. Substituting this in Eq. (j5.ip and reading the 0(e) term, we obtain 



d^i 9J"cff(fp,'C2) aj"eff(rp,/:^^ 

= — *3Tl 4- — 



dr 



eri 



e=0 



(5.3) 



where 5^ denotes a linear variation with respect to e (holding tq fixed). Since SeiC^) — for geodesies [see Eq. (|2.4[) ]. 
we obtain 



with 



d^n 

dT^ 



dr„ 



e=0 



M{ro - 6M) 
r3(ro-3M)' 



(5.4) 



(5.5) 



Equation (|5.4p tells us that the orbit is stable under small-e perturbations whenever > 0, and is unstable under 
such perturbations when < 0. The innermost stable circular orbit is identified by the condition = 0, giving 
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Hsco = 6M. We also find that, at 0(e), the radial motion is harmonic in r: Integrating Eq. (|5.4I) with the assumption 
of a periapsis passage at r = 0, we obtain ri — —tq coscj^t and hence 



rp(T) = ro(l — e cos w^t) + O(e^). 



(5.6) 



Next, we consider the 0(/i) correction to the orbit caused by the conservative piece of the SF. The equations of 
motion become 



— 1 TTicons 



dC 



— 1 TTicons 



(5.7) 



df2 



(5.8) 



where hereafter overtildes indicate quantities associated with the SF-corrected orbit (which is no longer a geodesic). 
Here we have defined the SF-corrected energy and angular momentum parameters £{f) and C(t) (in general no longer 
constants of the motion) through 



dt 



p 

df 



d(pp 
df 



f2 



(5.9) 



in analogy with Eqs. (|2.1[) . We assume that the orbit remains bound under the effect of the conservative SF, with 



< '^p('^) ^ '^max: ^ud ouce again define p and e as in Eq. 



, replacing r„ 



r, 



min and ^max 



(we leave 



e and p untilded for notational brevity). Without loss of generality we take rp(f = 0) = fmin- 

The radial component F^^^^ is an even, periodic function of r along the geodesic Xp(r) [recall Eq. (12.33^ ]. and 

hence also, at leading order in fi, an even, periodic function of f along the perturbed orbit ip(T). [This is because 

ip(r) — Xp{t) cx O(ijl) while FJons is already ©(/x^).] Since fp(T) too is even and periodic in f (and monotonically 

increasing between fmin and fmax), we may express F^^^^ as a function of fp only, for given p, e: -FJons — -^cons('^p;-P) s). 

In Eq. ()5.8|) the quantities fp(T), d^fp/df^ and -Fc'ons all periodic and even in r, and we conclude that C, too, is 

periodic and even. Hence, we may write C — C{rp]p, e). 

Let us now specialize to a slightly eccentric (SF-perturbed) orbit. Working through 0(e), we write rmin — ''o(l ~ e) 

and Tmax = ro{l + e), where ro[= pM + 0{e^)] is the radius of the circular orbit about which we perturb. For this 

orbit we write fp(T) = rg + efi(r) + 0{e'^) [as in Eq. (j5.2p ]. and we have Fc'ons = ^cons(^p; ^Oi e) and C = C(fp] tq, e). 

At 0(e°) (i.e., at the circular-orbit limit) C is constant along the orbit from symmetry, and solving Eq. (|5.8p with 

d^fp/df^ = immediately gives 



P - 



Mrl 
ro - 3M 



1 



— ^^0 



(5.10) 



where hereafter subscripts "0" denote circular-orbit values. In particular, we denote by Fq the circular-orbit value of 
(omitting the label 'cons' for brevity). Then, at 0(e), Eq. (|5.8p yields 



zpr 
cons 



d^fl 
'df2 



dfr, 



ef 1 -I- 



dFcsiro,C'^) 



fp=ro 



9(£2) 



"e^ cons' 



(5.11) 



where we have used 6efp — efi. To evaluate SeC^ and 5eF^^^^, we note that the two quantities depend on e both 
implicitly, thought rp(r;ro,e), and explicitly. However, as we showed in Ref. [s^, the explicit linear variation of 
and -FJons with respect to e (with fixed vq and Vp) vanishes at e = 0. Hence we may write S^C^ = efidC^ /dvp and 
JeFJons — efidF^^^g/dfp (where the derivatives are evaluated at e — 0), resulting in that Eq. (|5.1ip takes the form 



(5.12) 



with 



-2 



df„ 



J"cff (f p , ^ (f p ) ^ ) + M ^ -Picons (^P ) 



(5.13) 
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To obtain a more explicit expression for the shifted radial frequency w^, we next expand C and fJons ^ through 
0(e). First, Solving Eq. (|5.12p with the initial condition fp — fmin, we find fi = —VQCOsCorT. Then we expand 
Kons = + eri (rfFe'ons/drp),_^^,^ + 0(e2), giving 

= + ei^[ coscD.f + 0{e\ (5.14) 

where we have denoted Ff = — ro {dF^^^^/df-p)^ In a similar manner we expand £ = £o + efi{dC/dfp)f^=ro + 
O(e^), which, in conjunction with Eq. ()5.7p . gives 

F™"'' = et:^,^^ sintL.f + 0(6^), (5.15) 

where = iJ.ro{d£/ dfp)f_^=rg- Finally, using Eqs. (|5.14p and (|5.15p and substituting for Co from Eq. ()5.10p . Eq. (|5.13p 
becomes 

- ^2 _ 3(ro - 4_M) ^ _ 2 ^^^^^ _ 3M)a.-^FJ 



ro(ro-3Af) ro 



r3(ro-3M) 



3r§(ro - 4M) rg(ro - 3M) 2(ro - 3M)^M{ro - 3M) ' 



(5.16) 



This formula describes the 0(/i) conservative shift in the radial frequency off its geodesic value. Note that it requires 
knowledge of the SF through 0(e) (knowledge of the circular-orbit SF is not sufficient). 

The perturbed ISCO radius, r = fjsco, is now obtained from the condition w^(ro = Hsco) = 0. Namely, fjsco is the 
value of To that nullifies the expression in square brackets in the second line of Eq. (|5.16p through 0(/i). Note that 
in this expression we are allowed to substitute rg = risco ~ 6Af in all SF terms, since such terms are already 0(/i) [so 
the error introduced affects fisco only at 0(/i^)]. Thus, we readily obtain 

^^isco = ^isco QAI 

= {My^i)(216FS,,-m8F^,, + V3M-^F^,,) (5.17) 

through O(^), where we have denoted Fqj^ = Fq (t-q = 6M) and similarly for F[^^, F^^^. 

Since the coordinate ISCO shift Arisco is gauge dependent (just like the SF itself), it is not very useful as a 
benchmark for comparisons. Instead, we now consider the (SF-corrected) circular-orbit azimuthal frequency, 

n^^ = ^^, (5.18) 

dt dtp/df 

which, as discussed in jl9| . is invariant under all 0(/i) gauge transformations whose generators respect the helical 
symmetry of the circular-orbit configuration. Using Eqs. (15. 9p and (|5.10l) we obtain 



n = n 



ro(ro-3M) • 



(5.19) 



where O = {M/r^y/"^ is the geodesic (no SF) value. Evaluated at — fisco, the SF-induced frequency shift Afi = Q—U 
reads [through 0{fi)] 



1 



63/2 M 



, 27M 



4M 2n 



(5.20) 



Equations analogous to (|5.17p and (|5.20l) were obtained by Diaz-Rivera et al. [l3| in their study of scalar SF effects. 

Despite the fact that Afijsco is gauge invariant (in the above sense), care must be exercised in interpreting the 
quantity expressed in Eq. (j5.20p . As pointed out in p6j . the Lorenz-gauge metric perturbation has the somewhat 
peculiar feature that its tt component (in Schwarzschild coordinates) does not fall to zero as r — ^ oo, but instead 
htt — — 2a(= const), with a = /^[ro(ro — 3M)]~^/2^ This peculiarity can be removed simply by "rescaling" the time 
coordinate ast^i={l-\- 2aY^^t ^ (1 -|- a)t [neglecting terms of ©(/i^)]. The angular frequency Vt, whose definition 

has an explicit reference to t, will be modified under such rescaling as i7 $7 = (1 — As explained in Ref. 

[l^ . and recently emphasized by Damour [2l| . in order to compare with the frequency derived in a gauge in which 
htt admits the ordinary asymptotic fall off [such as the Regge- Wheeler gauge, or the gauge associated with the 0{im) 
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part of EOB theory], one must not use the "Lorenz-gauge" frequency O, but rather the i-rescaled frequency (l. For 
practical reasons, therefore, we also give here the "i-rescaled" version of Eq. (|5.20|) : 



1 



63/2 M 

where we have used a{ro = 6M) = (/i/M)/-\/T8. 



AM 



27M 
2(1 0- 



18M 



(5.21) 



B. Numerical method and results 



In view of Eqs. (|5.17p and (j5.20p . the task of calculating Arjsco and Ailisco amounts to obtaining numerical values 
for the three coefficients Fq-^, F^-^ and F^jg. The first coefBcient is easily obtained: Fq-^ is just the SF along a strictly 
circular geodesic with radius = 6M, and we already computed it in Paper I using our circular-orbit code (it is one 
of the values listed in Table V of [l^). Here we repeat this calculation at greater numerical precision, obtaining 

Fo'ig = 0.0244665(1) n/M^. (5.22) 

This is consistent with Berndtson's fl^ result of 0.024466497 /x/M^. 

The computation of F^^^ and F^-^^ is much more delicate, as it requires to resolve numerically the small variation in 
the SF under a small-e perturbation of a circular orbit [recall Eqs. (15.141) and (|5.15p ]. We are not helped by the fact 
that this variation need be evaluated at (p, e) — (6,0), which is a singular point in the p-e plane (see below). The 
subtlety of the computation task calls for an extra caution, so, as a safeguard measure, we devised and implemented 
two completely independent strategies for evaluating F^-^^ and -Ffjg. The first, more direct approach (method 7) involves 
an evaluation of the SF along a sequence of eccentric geodesies approaching the ISCO along a suitable curve in the 
p-e plane; the required SF coefficients are then extracted as certain orbital integrals, extrapolated to e — (see 
below). The second strategy {method II) involves an expansion of the field equations themselves about a circular 
orbit, through 0{e). In what follows we describe each of the two methods and their outcomes in turn. 



1. Method I: extrapolation in the p-e plane 
From Eqs. (|5.14p and (|5.15|) we obtain 

/"TT / LOj. 

F^i, = lim lim/'[(p,e), Fi''(p, e) = 2w,.(e7r)-M FJ^^^ cos w^rdr, (5.23) 

p^6 e^O Jq 



and 

F^, = hm hm F^(p,e), F^{p,e) = 2{e7r)-^ / F™- sin t^.r dr, (5.24) 

^ e^O ^ ^ Jo 

where w^, r, F^^^^ and F^°^^ are the values corresponding to a geodesic with parameters p, e (hence we were allowed 

to remove the tilde symbols off ojr and t). It may be noticed that, formally, the quantities F^{p, e) and F^{p, e) are 
inverse Fourier integrals describing the first w^-harmonic of FJ^j^g and F™"^, respectively. The latter two quantities, 
recall, are both periodic in r along the geodesic, with F^^^^ being even and F™"** being odd in t. 

The order of the limits in Eqs. (|5.23p and (I5.24p is very important: Since F^j^ and F^j^ are defined through an 
expansion about a stable circular orbit (e = 0), we must first take the limit e — >■ before taking p — >■ 6. In practice, 
however, it is more computationally economical to approach the point (p, e) = (6, 0) along a certain continuous curve 
in the p-e plane, rather than having to extrapolate to e — > along several different p=const lines and then extrapolate 
the resulting data again to p — > 6. In doing so, however, we must choose our curve carefully, since the limiting point 
(6,0) is known to be a singular one. [A simple manifestation of this singularity is the fact that the rate at which the 
radial frequency Ur vanishes at the limit (jo, e) — >■ (6, 0) depends upon the direction in the p-e plane from which this 
limit is taken — see, e.g., Eq. (2.36) of [Sg.] In particular, the curve must always "stay away" from the separatrix 
p — 6 + 2e, where the e-expansions (|5.14p and (|5.15p [on which Eqs. (|5.23l) and (|5.24p rely] are meaningless. As 
discussed by Cutler et al. in 3^ (in a slightly different context), in order for the e-expansion to hold, one must require 
not only e <C 1 but also e <C p — 6. For our purpose, the point (6, 0) must be approached keeping e ^ min{l,p — 6}. 

Here we pick the curve p = 6 + for which e <C 1 automatically implies e <C p — 6. We select a set of p 
values approaching p — ^ 6, and for each value we use our code to calculate the SF along an eccentric geodesic with 
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parameters {p, edp)), where edp) = (p— 6)^. For each of these parameter-space points we then construct the quantities 
F[{p, ec{p)) and F^{p, edp)) defined in Eqs. (|5.23p and (|5.24l) . by numericahy integrating the SF data along the orbit. 
The results are shown in Fig. ^ Finally, we obtain the desired coefficients F^-^^ and F^-^^ by extrapolating the numerical 
values of F[ and .F^ to p = 6. We do this, in practice, by fitting a cubic polynomial to the numerical data, writing (for 

example) F[ = uq + ai(j) — 6) + ■ ■ ■ + a3{p — 6)^ and taking the value of the fitting coefficient oq as our approximation 
for F^jg. The results are 

F[.^ = 0.0620(5) Ai/M^, 

F^-^ = -1.066(1)/! {method T). (5.25) 

Here, as elsewhere in this work, a parenthetical figure indicates the uncertainty in the last displayed decimal place 
due to numerical error [so, e.g., 0.0620(5) stands for 0.0620 ± 0.0005]. 

The error in the above calculation is estimated as follows. For each value of e considered, we first estimate the 
numerical error in each of Ff and F^ as the difference between the value obtained with the finest numerical grid used 
and the value obtained with a coarser grid (of four times the cell area). We calculate this difference for each of the 
Z-modc contributions to Fl and F^, and conservatively take the total (Z-summed) error as the sum over the moduli of 
the individual l-moAe errors (adding to this the large-Z tail extrapolation error). The "error bars" thus obtained are 
those displayed in Fig. [HI In the final step we fit a cubic polynomial curve to each of the Ff , F^ data sets, using the 
above error bars as a fitting weight. We take the standard fitting error in the constant term (the above coefficient oq) 
as our estimate for the error in Fy^^.F^-^^. It is this estimated error that we indicate in Eqs. (j5.25p . 
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FIG. 9: Derivation of the SF coefficients Fy^^ and F^^^ entering the ISCO-shift formula (|5.17|l . using method I. Plotted are the 
numerical values of F[{p, e) and F^{p, e) [see Eqs. 1)5. 23|) and (|5.24|) ] for a sequence of points in the p-e plane approaching the 
ISCO along the curve p — 6 + ^/e. Error bars indicate the estimated numerical error, which is evaluated as explained in the 
text. The solid curves are weighted polynomial fits used to extrapolate the values of F[ and F^ to p ^ 6. The extrapolated 
values are the desired coefficients Fi^^ and F^ig. 



2. Metfiod II: e-expansion of ttie field equations 

Our second procedure for evaluating F^jj, and F^^^ is based on a systematic expansion of the field equations (and 
the SF) through 0(e), and is similar, in principle, to the method used in Ref. [l3| for the scalar- field SF. The idea is 
to consider a slightly eccentric geodesic, with rp(T) = ro(l — ecosuirr) where e ^ 1 and ro=const [as in Eq. (|5.6|) ]. 
and for this geodesic expand the source term in the perturbation equations (j2.20l) through 0(e) as 



{t, r)S{r - rp) = {t, ro)S{r - ro) + e S['\t, rQ)S{r - ro) + T^\t, ra)5'{r - tq) 



(5.26) 



Here we have omitted the indices Im for brevity and used a prime to denote d/dr. The various expansion coefficients 
are given by 



(,) _ a5«(i,rp) 



de 



e=0 



cii) 

= Jq ro cos UtrT. 



(5.27) 



e=0 
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FIG. 10: Derivation of Fy^^ and F^is, using method II. Plotted are the numerical values of i^c'os(p) a-nd F^'"(p) [see Eqs. (|5.33|) 
and (|5.35p ] for a sequence of points approaching the ISCO along the p-axis in the p-e plane. Error bars indicate the estimated 
numerical error, which is evaluated as explained in the text. The solid curves are weighted polynomial fits used to extrapolate 

to p — > 6. The extrapolated values are the desired coefficients F[ig and F^^^. 



The coefBcients S^ ^ (for i — 1, . . . , 10) are evaluated by substituting the above expression for rp^r) in Eqs. (|Blip - 
(iB2?l . along with u'' = erow^ smujrT+0{e^), £ = {ro-2M)[ro{ro-3M)]-^/^+0{e^) and£ = rQ{ro/M-3)-'^/'^+0{e'^), 
and then expanding through 0(e). Explicit expressions for the S^"^ 's are given in Appendix|Fl We also formally expand 
the metric perturbation functions /i^*)'™ in the form 

U') (t, r) = h''^ {t, r) + eh['^ {t, r) + 0{e^) (5.28) 

(again omitting the indices Im for brevity), and consequently single out the 0(e°) and O(e^) pieces of the perturbation 
equations (|2.20p as 

ahl;^+M%hi'^ = S^'^6{r~ro), (5.29) 
Dh^;^ + M%h['^ = Si'^6ir - ro) + r|^'<5'(r - ro). (5.30) 

Notice that the source terms in these equations are evaluated along a circular geodesic (of radius ro). The function 
hg"^ is just the physical Lorenz-gauge perturbation from this circular geodesic. The function ft,^'-* is not a physical 
perturbation; in particular, it is discontinuous across the orbit (due to the S' term in its source). 

For our purpose we need to solve Eqs. ((5?29)) and (lOO)) for tq = 6M. While there is no problem solving for Hq^ 
(as we in fact already did in Paper I), solving for ft,^'-* is numerically subtle, since the source functions become 
singular as cj^ — J' [cf. Eqs. (IF2p ~ (|F12p in Appendix IF] . Hence, as in method I, we use extrapolation: we first solve 
Eq. (|5.30p for a sequence of tq values approaching 6M, and then extrapolate our data to ro — 6M. 

A second difficulty arises because our TD algorithm assumes that the numerical variables are continuous on the 
particle's orbit. This is no longer the case for the functions h^i \ which suffer finite jump discontinuities across the 
orbit. This requires an adaptation of the junction conditions implemented in the finite-difference scheme. We describe 
the suitably modified junction conditions in AppendixlF] We use our amended TD algorithm to solve for the functions 
/ij*'' via time evolution. 

Once the functions h^^ {t, r) and /i^*' (i, r) are at hand we can proceed to construct the SF along the slightly eccentric 
orbit through 0(e). Here we need to be cautious: When we evaluate /i'^*'(i,r) (and its derivatives) along the slightly 
eccentric orbit rpir), we must take proper account of the 0(e) contribution coming from the term /ig*^(t,rp) via the 
expansion of Vp in e. Recalling Eq. (j5.28p we have, in fact. 



U'\t,rp) = ^;\t,ro)+e[hf\t,ro)-hl;lit,ro)roCOSUJrT^ +0(6^), (5.31) 
with similar expansions applying to the r and t derivatives of the perturbation along t'p(t): 

/i«(t,rp) = h^^^,{t,ro) + e \h^^l{t,ro) - 4%{t,ro)ro cosuJrr] + 0{e^), (5.32) 
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where a = r or t. Recall that the derivatives of Hq"^ and h^i \ as well as the function h^^ itself, are discontinuous at 
rg and hence defined only in the sense of a directional limit r — > (however, for simplicity we choose here not to 
reflect this with a suitable notation). Notice also that the first derivatives of the metric perturbation h^"^^ along the 
slightly eccentric geodesic involve second derivatives of the function /iq*-* along the circular geodesic (an rr derivative 
for h^r and an rt derivative for ft. *'(■'). The necessary data for our calculation therefore includes the functions /ig*-* and 

themselves, as well as the derivatives Hq^^, Jiq^^ hi^r^ ^I't, ftg^r-r ^-^id ^o^rtJ evaluated (in a one-sided fashion) at 
the circular orbit with radius rg = 6M. 

To construct the SF through 0(e), we substitute Eq. (|5.28l) in Eq. (I2.25p . and evaluate the outcome along the slightly 
eccentric orbit ?'p(t). We remind that the functions J- appearing in Eq. (I2.25|) are given explicitly in Appendix [Q in 
terms of the perturbation h!^^\t,r) and its r and t derivatives. We next expand Eq. (|2.25p in e through 0{e) using 
the expansions (|5.31|) and (|5.32|) and the above e-expansions for , E and C. Keeping only the 0{e) terms, we 
proceed through the mode-sum regularization procedure as described in Sec. Ill El and construct the 0(e) pieces of 
the conservative SF components FJ^^g and F^^^^. We find that these 0{e) pieces contain only terms proportional to 
either coswr''" or sinwrT. More precisely, we obtain 

0(e) piece of FJ'o,,, = ei^J^^ cos w^r, (5.33) 
0(e) piece of i^™"" = ew^i^f sin w^r, (5.34) 

where the coefficients F^^^ and F^™ are constructed from hi\ hi^a-, h^^^ and /lo'.ro-- (The exphcit form of these 
coefficients is rather complicated and will not be shown here. The coefficients can be readily evaluated using computer 
algebra tools.) The 0(e) part of the orthogonality condition u°'Fa = then also gives 

0(e) piece of F™"" = ew^F^'" sincj^r, (5.35) 
where F^™ is related to the numerically-computed coefficients F^™ and Fq through 

--T^ i^oFr + roFS) . (5.36) 

Comparing Eqs. (jOH)) and with Eqs. (jgUl) and (I5J51) . one finally identifies Ff = F^^^ and F^ = F^'". 

We have implemented the above procedure for a dozen or so radii rg near 6M, and for each radius obtained the 
values of F^ and F^. We then extrapolated these values to rg = 6M by fitting a cubic polynomial (see Fig. [TU]). We 
obtained 

F^-^ = 0.062095(1) At/M^, 

F^i, = -1.0665(8) {method IT). (5.37) 

These values are in agreement, within the estimated numerical accuracy, with the values obtained using method I 
[compare with Eqs. (|5.25l) ]. Using method II we were able to explore orbits much nearer to the ISCO, resulting in 
a much improved accuracy. Of course, the agreement between our two independent calculations provides significant 
reassurance. 

With the SF coefficient values given in Eqs. and ([07)) . formulas (|CT7)) . ([g?^ and ([OT|) finally yield 

Ari,eo = -3.269(2)^, 

= rJisco X 0.4869(4) ^/M, 
AJ7i,eo = riisco X 0.2512(4) A*/M, (5.38) 

where riisco = (6^^^ Af)~^ is the unperturbed (geodesic) value of J7 at r = 6Af. The values obtained here are slightly 
more accurate than the ones we give in our Letter 33] [Arisco = —3.269(3) /i and Afiisco = ^^isco x 0.4870(6) fi/M], 
an improvement made possible by the inclusion of additional numerical data points in our analysis. 



VI. CONCLUDING REMARKS AND FUTURE APPLICATIONS 



This work marks a new frontline in the program to model realistic two-body inspirals in the extreme mass-ratio 
regime. For the first time we are able to calculate the full [0{fi^)] gravitational SF across (essentially) the entire 
parameter space of strong-field bound geodesies in Schwarzschild spacetime. This work also represents a first complete 
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end-to-end implementation of a range of computational techniques which were developed gradually over the past 
decade: mode-sum regularization scheme [l^, the 1-l-lD Lorenz-gauge perturbation formalism [2f|, and the method 
of extended homogeneous solutions [3l|. As the reader may appreciate, the underlying computational challenge is 
rather daunting, given the complexity of the field equations, the technical subtleties involved in dealing with the 
delta-function source, the high computational cost, and the need to patch together different techniques in both the 
time and frequency domains. Our eventual working code is of considerable complexity and took over two years to 
develop and test. 

Following is a summary of the various tests which helped us establish confidence in our code's performance, (i) 
The mode-sum regularization procedure is self-validating, in the sense that it is extremely sensitive to errors in the 
computation of the perturbation multipoles (especially the high-Z ones, which are most computationally demanding). 
If the regularized mode sum shows the expected fall-off behavior at large I, this by itself is a strong indication that 
the high-Z modes were calculated correctly, (ii) The code reproduces the known results in the circular-orbit case; 
these results are now confirmed by 3 independent analyses [ill ■ (iii) Our evolution code reproduces the correct 

asymptotic fluxes of gravitational-wave energy and angular momentum, as verified by comparing with results in the 
literature, (iv) The work done by the dissipative piece of the computed SF is found to precisely balance these fluxes, 
(v) The value of the ISCO frequency shift derived from our SF seems consistent with the value derived in FOB at 
3rd post-Newtonian (PN) order: Damour recently showed that the latter is about 72.5% of the SF value, with the 
difference likely attributed to higher-order PN terms 21]. 

In principle, our code can return the SF along any bound geodesic in Schwarzschild geometry, although in practice 
computational cost may becomes prohibitive when the orbital period is too large (i.e., for very large p and/or e close 
to unity). The "workable domain" of our code using a current-day standard single-processor desktop computer is 
roughly < e < 0.5 and p < 20M if a fractional accuracy of < 10~^ in the SF is sought. The current algorithm 
incorporates an explicit reference to the radial frequency parameter, so it cannot be used to tackle unbound orbits. 
However, it may be adapted with moderate effort to handle unbound orbits (including orbits below the last stable 
orbit) as well. 

Even within the above "workable domain" , the current code is discouragingly slow. It takes a few days to compute 
the SF along a single strong-field geodesic, which makes it impractical to cover the entire parameter space of inspirals at 
sufficient resolution (having in mind the development of theoretical gravitational waveform templates for astrophysical 
inspirals). There are, however, various ways in which one may improve the computational efficiency and speed. 
Most obvious, one can use distributed computing — our algorithm should be easily amenable for distribution on a 
cluster, since different I modes can be calculated in parallel. Other natural approaches include the use of mesh 
refinement (see Thornburg's recent report [H^l) and/or higher-order finite-difference schemes. One may also seek to 
reduce the computational cost attached to the initial stage of the numerical evolution (when spurious initial waves 
dominate) by iteratively improving the initial conditions for the evolution — this idea is already being implemented 
successfully in a 2-l-lD framework [551] . Other ideas represent a more significant deviation from our approach: (i) 
work entirely in the frequency domain, making full use of the method of extended homogeneous solutions [31| (this is 
likely to prove most efficient with smaller eccentricities); or (ii) abandon finite differencing altogether and instead use 
finite elements or other pseudospectral techniques (sgI . \5T\ , benefiting from their natural flexibility in accommodating 
multiple lengthscales. 

Nonetheless, some interesting applications are already possible with the current version of our code, and we have 
presented one of them in Sec. |Vl Our computation of the ISCO frequency shift represents the flrst physically- 
meaningful new result coming out of the SF program, and it has already informed both Numerical-Relativistic 
calculations [13, HH and EOB/PN theory [2l[. Further ideas for SF/EOB synergy were recently discussed by Damour 



21[. In particular, Damour showed that a computation of the (gauge invariant) conservative SF correction to 



the precession rate of the periapsis, for slightly eccentric orbits, will give access to the presently unknown 4PN (and 
possibly higher-order) parameters of FOB theory. We are currently working to extract the necessary SF data to 
facilitate this calculation [these are essentially the coefficients Fq, F[ and of Eqs. (I5.14[) and (I5.15[) . as functions 
of ro]. Our preliminary results show excellent agreement with the analytic FOB predictions at 2PN and 3PN, and we 
are hoping to publish these findings elsewhere |58|. 

More generally, our code allows to tackle the calculation of post-geodesic [0(/i)] precession effects at any eccentricity, 
not necessarily small. This would not only provide a handle on the 'Q' function of FOB theory, but can also directly 
inform the computation of conservative effects in inspiral trajectories. Information on the periastron advance as a 



function of p, e could, for example, be incorporated into an orbital evolution scheme a la Pound & Poisson [2j]. We 
are currently investigating this direction. 

Finally, we comment briefly on the possibility of extending this work to the Kerr case (a more elaborate discussion 
of possible strategies for attacking the Kerr problem can be found in [^). The framework of the current code, i.e., 
numerical evolution in 1-l-lD is not directly applicable in Kerr spacetime, because the perturbation equations in Kerr 
cannot be separated into harmonics in the time domain. It is possible to tackle the field equations in 2-l-lD (i.e, 
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separating the perturbation into azimuthal m-modes only) or in full 3+lD, and there has been considerable progress 
in that direction in the past two years — although work so far has been restricted to the toy problem of a scalar field 
in Schwarzschild geometry [111, [sl, Isol l60l| . Schemes for regularizing the SF directly in 2+lD or 3+lD have been 
proposed and recently implemented jSSl l60l. Isij . Alternatively, one may attempt to tackle the field equations in 1+lD 
by properly accounting for the coupling between different Z-harmonics, and a similar strategy may be applicable in 
the frequency domain too. A first calculation of the SF in the Kerr case (using the frequency-domain approach) — for 
a scalar charge in a circular equatorial orbit — will be presented in a forthcoming paper [63 |. 
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and for permitting us to include these data here. We are grateful to Thibault Damour for useful comments on the 
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Sasaki and the Yukawa Institute for Theoretical Physics in Kyoto (where part of this work took place) for their 
hospitality. NS acknowledges support from Monbukagaku-sho Grant-in-Aid for the global COE program "The Next 
Generation of Physics, Spun from Universality and Emergence" . 



We prescribe here the construction of the various components of the (trace-reversed) Lorenz-gauge metric per- 
turbation hai3 in terms of the 10 time-radial scalar-like functions introduced in Eq. (j2.19p . In the following 
yZm _ Y''"^{9,(p) are the standard spherical harmonics, / = 1 — 2M/r, and for brevity we suppress the multipolar 
indices l,m in /j^')'™. Recall that the functions — our basic numerical evolution variables — are obtained by 
solving the coupled set (|2.20p . 

The metric perturbation is reconstructed through 
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Appendix A: Algebraic reconstruction of the metric components in the Lorenz gauge 
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The angular functions appearing in these relations are defined as 
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We note that for Z = 0, 1 we have h^'^'^^^ = identically, and that for / = we have additionally /i(4,5,8,9) _ q_ 
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Appendix B: Field equations and gauge conditions for the perturbation functions h''''^'"^ (r,t) 



We give here explicit expressions for the various terms appearing in our 1+lD field equation (|2.20l) . In what follows 
/ = 1 — 2M/r, /' = 2M/r^, is the standard tortoise radial coordinate defined through dr/dr^ — f{r), and v = f+r*. 
We also denote A = (/ + 2)(Z - 1). 

The terms in Eq. ({2:201) read 
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(BIO) 



These expression are the sarae as those given in Appendix A of paper I, although we write them here in a shghtly 
different form, more amenable to discretization in v, u coordinates. 
The various source terms 

g(i)lm ((2:201) are given (referring to Sec. Ill A] for notation) by 
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The 1+lD field equations (|2.20p are supplemented by 4 elliptic "constraints" stemming from the Lorenz-gauge 
conditions (|2.15p . These constraint equations read 
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Appendix C: Construction of the full-force spherical-harmonic modes 



Following are the explicit values of the functions Tf'^V^ appearing in Eq. (|2.25p : 
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The various functions ^re those appearing in Eq. ([2?24)) . and below we give these functions explicitly for a — t^r 
(the values for are not needed in this work). The values of the various Z, m-dependent coefhcients a, /3, 7, e, C, C in 
the above expressions are given in Eqs. (|D3p - (|D9P of Appendix IdI In what follows we use the notation C = C/r-p, and 
for brevity we suppress the indices Z,to as well as the subscript ±; it is to be understood that the t and r derivatives 
in the following expressions are taken from either "outside" or "inside" the orbit, yielding, in general, two different 
one-sided values which correspond to the suppressed + or — subscripts. 
For the t component we have — 
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For the r component we have 
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Appendix D: Useful identities 

The following identities are used in deriving Eq. p.25p for the full-force modes. In these relations F'™ are the 
standard spherical harmonics, and the identities are valid for any values of /, m. 
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Here the various coefficients are all constructed from 
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Appendix E: Jump conditions for the perturbation modes and their derivatives 

As explained in Sec. IIII Bl our finite-difference scheme makes use of formal jump conditions for the perturbation 
modes and their (first through fourth) derivatives across the particle's orbit. In this appendix we derive the 

necessary conditions. Our derivation refers to a specific (yet generic) point xq along the orbit, with known coordinates 
(ro,to) or {uo,vq), and velocity components uq = dup/drj^^ and wo = dvp/dT\^^, where t is proper time along the 
orbit. We will use the notation = A{xq) — A{xq), where A{x^) are the values of a 1+lD field A(r,t) calculated 
by taking the limits t ^ to and r — >■ r^. In this appendix we shall omit multipolar indices Im for brevity. 



1. Continuity condition for /i''^ 



The Lorenz-gauge perturbation modes /i^*^ are all continuous at any point along the orbit. This can be verified, 



for example, by noticing that the distributional form /i'*^ = h''^^ {r,t)9[r — rp{t)] 

solution of the perturbation equations (j2.20p only if the homogeneous solutions hj^ 
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2. Jump conditions for the 1st derivatives 

Let us we re-express the field equations (I2.20p in the form 

/.W +pW = J dTS^^HMr))S{u-u^{r))S{v-v^{T)), 

where 
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with A^^(]) and S'*^*^ being the quantities given in Appendix iBl and V{r) being defined in Eq. (j2.2ip . Now consider 
formally integrating Eq. (|E2p along the ray v = vq over an infinitesimal interval {uq — e, uo + e) across xq. The integral 
of the left-hand side yields 



as e — 0, since P^^> is bounded. Similarly integrating the right-hand side of Eq. (|E2[) one finds 



The desired jump condition is therefore 
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Similarly integrating the field equations along u — uq over an infinitesimal interval (wq — e, -I- e) one also obtains 
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3. Jump conditions for the 2nd derivatives 

It is straightforward to derive the jump condition in the mixed ww-derivative. From Eq. (jE2[) we immediately have 
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as the source term is supported only on the worldline. Recall P'*^ involves the perturbation /i^'^ and its first derivatives 
only, and so the jump in P^'^ can be readily calculated from the jump conditions obtained above. 

To derive a jump condition for the w-derivative, we first take the v derivatives of the field equations (|E2[) . and 
then integrate with respect to u along v = vq over {uq — e,uo + e). We thereby obtain 



~ lim / du 

e^O 



uo~e 



Now, for P^*) we may write 



P 



(E9) 



(ElO) 



in which pf^ are smooth and where Wp [u) is the value of v at which the outgoing ray of retarded time u intersects 
the worldline. Using this form we obtain 
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as e — 0, since P±^\ are bounded. Lastly, for the integral involving S'^^'^ in Eq. (|E9[) we have 

du / dr5«(xp(T))(5'[i>o - VpirMu - Up(r)] 
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where tq is the value of r at xq, and 6± are the values of r at which the two outgoing rays of constant retarded times 
Wo ± e cross the worldline. Combining the above results, we arrive at 
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The jump condition for the uu-derivatives is derived in a similar fashion. The result is 
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4. Jump conditions for the 3rd derivatives 

The jumps in the two mixed third derivatives are readily obtained by considering the derivatives of the field equations 
E2p with respect to v and with respect to u: 
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Here the jumps in Pli^ and P.i'^ can be obtained from the jumps in the first and second derivatives of /i*-*^ derived 
above. 

To obtain the jump condition for the vvv derivative, we differentiate the field equations (jE2l) twice with respect to 
V, and then integrate with respect to u along v = vq over {uq — e,UQ + e). This yields 



= lim 

e->0 



uo+e 



du 



P^liu,va) + / dTS^'\xp)S"[vQ - Vp{T)]5[u - Up(r)] 



(E17) 



Using Eq. (jElOp again, we obtain 



uo+e 



P^l{u,vo)du = I du{2(p^l{u,vo)~P%{u,vo))5[vo-Vp{u)] 

+ [p^\u,vo) - P^'Hu,vo)) S'[vo - v^{u)] 

vaWivo - Vpiu)] + vo)9[vp{u) - vq]^ 



p{i) 



(E18) 
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as e — >■ 0. Here we have introduced v'^ = dvp{u)/du and v'q = v'p{uo) = vq/uo] and similarly, defining Up{v) as the 
value of u at which the incoming ray with advanced time v intersects the worldline, we introduced u'^ = dup(v)/dv 
and Up = d"^ u-p{v) I dv^ , with Uq = Up(wo) = uq/vq and Uq = Up(wo) = ('"o'^o ~ uoVo)/vq. To evaluate the limit in Eq. 
(jElSP we have used 



uo+e 



du (P^;l{u. «o) - P%{u, «o)) 5[vo - Vp{u)] ^ {v'^) ' 



(as e — >■ 0), along with 



tio+e 



du (Pf (u,t;o) - P^"\u,vo)) S'[vo - vp{u)] 

du<[p'i^\u,v„)-PL'\u,vo)) K)-'-^ [{v'pr^5{u-up{v))] 



uo-e 



du <i 4- 

du 



p(0 

p(0 



p(i 



P 



The term involving S'^^^ in Eq. (jE17l) gives, upon twice integrating by parts, 

''P dT TP * 

as e — 0. Finally, substituting from Eqs. (jE18|) and (|E2ip in Eq. (|E17P we arrive at 

+W)-'[p'.;i+.sM„-«„-ik'^(v'^") 



du / dT5«(xp)<5"(z;o-«p(r))5(u-Up(T)) ^ 



The jumps 



p(i) 



Xq 



P 



and 



p(«) 



(E19) 



(E20) 



(E21) 



(E22) 



can be obtained in a straightforward way from the jumps in the first and second 



derivatives of the perturbation, obtained in previous steps. 

In a precisely analogous manner, differentiating the field equations (jE2p twice with respect to u and then integrating 
with respect to v along u — uq over (vq — e, wq + e) gives 



- 2K)- 



p('i 



p(0 



(E23) 



where = v''{uo) = (iiouo - woiio)/uo. 



5. Jump conditions for the 4th derivatives 

The jumps in the five 4th-order partial derivatives of /i^*^ are obtained in a similar manner. The results are 
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Here u'^' = Wp"(ro) = (uovl-'iu(iVQVQ + 'iuovl-uov\vo)/vl and v'^' = v'^'{ro) = {vouI-3voUoUq+3vouI-voUoUo)/u^ 
The jumps in the various P^'-* terms are related to the jumps in the 1st, 2nd, and 3rd derivatives of the perturbation, 
which were obtained in previous steps. 



Appendix F: Source terms and jump conditions for slightly eccentric orbits 

1. Source terms 

We give here exphcit expressions for the 0(e) source coefficients 5^;*^ defined in Eq. (|5.27p . These are needed in our 
discussion of the e-expansion method in Sec. |Vl We use the notation 



£q = £{p = ro,e = 0) 



ro - 2M 
v/ro(ro-3M)' 



Co ^ C{p = ro,e = 0) 



Mro 



(Fl) 



along with /o = 1 — 2M/ro, and omit multipolar indices Z, m for brevity. 

The coefficients S^'' are obtained by formally expanding the source terms 5*^*^ given in Eqs. (jBliP " (jB20p in powers 
of e through 0(e), assuming the e-expansion forms of the trajectory, 

rp(T) — ro{l — e cosw^r) + O(e^), fp{T) ~ (-^JipT H —e sincj^T + O(e^), 



where we wrote = Cq/tq. S^^ in then obtained as the linear variation of S^^'^ with respect to e [recall Eq. (|5.27p ] 
The result can be expressed in the form 



•^1 — 2^ '-'l.Tl 



n=±l 



where tOnm = nuJr + toWm and the various coefficients are 



c(l) 
c(2) 



(3) 



SI 



Cj(4) 



27r(ro - 2M)3 
for-^(rn-3M) 

27rfo ' 



ro - 6M + 2mn{ro - 2M) 



ro - AM + 2mn{ro - 2M)-^ 



SnimCoiro - 2M) 



ro - 4Af + mn{ro ~ 2M)^ 



ST. 



(5) AwmnLUrCoiro — 2M)'' 



£or$ 



(6) _ 27rCliro-2M) 

'-'l,n — 



c(7) 
o(8) 

o(9) 
'->l,n 

c(10) 



So4 

[lil + l)-2m']s[% 
SirCoiro - 2M) 



3ro - lOM + 2mn{ro - 2M)- 



ro - 4M + mn{ro - 2M)-^ 



4:TrinuJrCo{ro — 2My 



A'Ki'mCl{ro - 2MY 



YLA^), 



3ro - lOAf + 2mn(ro ~ 2M)-^ 



YLom, 



(F2) 

(F3) 
(F4) 
(F5) 

(F6) 

(F7) 

(F8) 
(F9) 
(FIG) 

(Fll) 

(F12) 



Here uJnm = ncur + mLL>^, with Yj^{0) and Yj^ g{0) denoting the spherical harmonics and their ^-derivatives evaluated 
a.t 9 = tt/2 and (p = 0. 

Notice that if we also write 



with % 



(F13) 



n=±l 
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then we may split the 0(e) part of the metric perturbation into two harmonic components, fi'f' — h}'i'^i + h\^'_i, each 
of which satisfying the field equation 



sfi5{r-ro)+Tri5'{r-ro) 



(Or// 



{n = ±1) 



(F14) 



Evidently, the source for ft,^*^ is simple harmonic in r, with frequency 



It is clear that the solutions ^xnh'-i^^ 



(and their derivatives) will inherit this simple harmonic dependence on r along the orbit. Hence we may write, for 
example 



where the factor /o£q ^ is simply [dtp/dT) ^. 



rp(T) 



— — *W„„j/oto ""l.n 



rp(T) 



(F15) 



2. Jump conditions 

Here we explain the derivation of the necessary jump conditions for the variables h^i^ associated with the 0(e) 
piece of the metric perturbation [recall Eq. (|5.28p ]. Unlike the full physical metric perturbation functions ft.*-*-' [or the 

0(e°) variables /ig*'], the functions lif' are in general discontinuous across the particle's worldline, because the field 
equation which defines them, Eq. (|5.30l) , is sourced by a derivative of a delta function. For our TD algorithm we need 

to derive analytic jump conditions for h\ itself as well as for its first to fourth derivatives along the orbit. Our task 
here will be somewhat simpler than in the case of the full perturbation /i^'^ (analyzed in Appendix |E]) thanks to the 

simple harmonic dependence of Ji'i^ on t along the orbit, expressed in Eq. (jFlSp . 
We start by reexpressing the field equations (jF14p in the form 



{dl - 9,2) h^^ ~ AP[]l ^ -4 5l:iJ(r - ro) + T^^S' {r - ro) 



where the definition of p[''\ is similar to that of P^^^ in Eq. (|E3[) . simply replacing /i*^*) — > 'li'i- Substituting 



(n = ±l), 



(■0 



(F16) 



(F17) 



and comparing the 5' {r — tq) and 8(r — ro)-terms on both sides of the resulting equation, we readily obtain the jump 
formulas 



Jo 



fo 



-E 



5,n 



4 
H 



jQ'l,n 



where Pg^j^ are the terms oc S{r — vq) in which are given by 



(0 



(F18) 

(F20) 

(F21) 

(F22) 

(F23) 
(F24) 

We may now easily obtain jump conditions for the second and higher r,. derivatives in a recursive manner: Formally 
differentiating Eq. (jF16|) k times with respect to r» we get the jump relations 
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(fc>0), 
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where the second equahty is due to the harmonic dependence of h^l\ on t along the orbit, expressed in Eq. (jF15l) . 
The jumps in the t derivatives are also obtained in a straightforward manner, by writing 



(F26) 



The desired jump conditions in h\ and their derivatives are finally obtained by adding up the jumps in the two 
n — 1 harmonics: 



(F27) 



n=±l 
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